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Local  scour  around  bridge  piers  is  a  significant  problem  in  the  design  and 
maintenance  of  bridges.  Knowledge  of  the  rate  at  which  scour  occurs  can  help  engineers 
in  the  areas  of  design,  monitoring,  and  scheduling  of  repairs.  The  temporal  variation  of 
clear  water  local  scour  around  a  single  circular  cylinder  is  examined.  Physical  experiments 
are  conducted  using  erodible  bed  and  fixed  bed  models  of  a  single  scour  event.  The 
erodible  bed  experiments  include  a  near  equilibrium  scour  experiment  and  two 
experiments  that  are  stopped  at  intermediate  times  during  the  scouring  process.  The  fixed 
bed  experiments  model  the  holes  created  during  the  erodible  bed  experiments  as  well  as  an 
initial  flat  bed  condition.  The  velocities  of  the  flow  around  the  cylinder  for  these  four 
cases  are  measured  using  an  Acoustic  Doppler  Velocimeter  (ADV).  Numerical 
simulations  of  the  flow  through  the  scour  holes  are  conducted  using  the  commercially 
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available  three  dimensional  Navier-Stokes  equation  solver  FL0W3D®.  The  resuhs  of  the 
simulations  and  the  measurements  from  the  physical  experiments  are  compared  to  justify 
the  use  of  the  software  for  further  numerical  experiments.  The  software  is  used  to 
produce  numerous  simulations  of  the  scour  event  at  intermediate  depths.  These 
simulations  are  used  to  develop  a  method  for  predicting  the  depth  of  scour  time  history. 
This  method  is  based  on  the  competing  mechanisms  of  the  erosive  and  resistive  actions  of 
the  hydrodynamics  and  the  scour  hole  geometry.  The  method  relates  the  volumetric  rate 
of  sediment  leaving  the  hole  to  the  kinetic  energy  at  some  fixed  distance  above  the  surface 
of  the  scour  hole  and  the  potential  energy  associated  with  removing  sediment  from  the 
bottom  of  the  scour  hole.  This  method  is  validated  through  the  analysis  of  a  second  scour 
event  that  has  entirely  different  flow,  sediment,  and  structure  properties. 
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CHAPTER  1 
INTRODUCTION 

1.1  Motivation 

On  April  5, 1987,  ten  people  were  killed  by  the  collapse  of  the  New  York  State 
Thruway  bridge  across  the  Schoharie  Creek.  On  April  1,  1989,  an  80  foot  section  of  the 
U.S.  51  bridge  over  the  Hatchie  River  near  Covington,  Tennessee,  toppled  into  the  river 
after  two  concrete  pile  bents  supporting  three  concrete  spans  failed.  Eight  casualties  were 
reported.  On  May  26,  1989,  two  spans  of  a  temporary  bridge  crossing  the  Great  Miami 
River  near  Miamitown,  Ohio,  collapsed.  Two  deaths  resulted  from  the  bridge  failure.  The 
cause  of  the  failure  of  each  of  these  bridges  was  identified  as  scour  around  the  foundation 
piles.  These  unfortunate  events  have  thrust  national  attention  on  the  pier  scour  problem. 
As  a  consequence,  there  is  a  strong  desire  world-wide  to  improve  scour  depth  prediction 
methods  for  the  design  of  pier  foundations  and  to  improve  methods  for  determining  the 
safety  of  existing  bridges  once  a  scour  hole  has  been  formed. 

The  problem  of  scour  around  bridge  piers  is  significant.  In  a  comprehensive  study 
of  bridge  failures  in  the  United  States,  Brice  and  Blodgett  (1978),  under  the  auspices  of 
the  Federal  Highway  Administration  (FHWA),  studied  383  bridge  failures  caused  by 
catastrophic  floods.  Their  findings  indicate  that  approximately  half  of  these  failures  were 
caused  by  local  scour.  Of  these  failures,  some  of  the  local  scour  was  attributed  to  ice  or 
debris  being  caught  on  the  piles  thus  increasing  the  local  and  contraction  scour.  A  large 
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portion  of  these  failures  resulted  from  erroneous  prediction  of  local  scour  during  design 
events.  The  report  states  that  the  damage  to  bridges  from  floods  between  the  years  of 
1964  to  1972  amounted  to  approximately  $100,000,000  per  failure. 

After  the  catastrophic  flood  of  the  Mississippi  River  in  1993,  Kamojjala  et  al. 
(1995)  published  a  report  detailing  the  damage  and  resulting  costs  of  bridge  failures 
associated  with  the  event.  Of  the  28  bridge  failures  that  occurred  on  Federal  Aid  Routes, 
two  of  the  failures  were  directly  attributed  to  local  scour  totaling  almost  $1,000,000  in 
damages.  Many  more  of  the  2500  bridges  requiring  FEMA  assistance  also  experienced 
damage  due  to  local  scour. 

These  studies  underscore  the  importance  of  prediction  of  local  scour  in  terms  of 
the  costs  both  monetarily  and  in  human  life.  It  is  obvious  that  better  understanding  of  all 
facets  of  the  local  scour  phenomena  is  of  vital  interest  to  the  engineer  responsible  for  the 
prediction  of  local  scour  in  the  design  of  pier  foundations.  The  FHWA  has  sponsored  a 
number  of  investigations  aimed  at  improving  methods  for  prediction  of  local  scour,  and  is 
also  involved  in  a  number  of  monitoring  programs  whose  goal  is  to  identify  problem  areas 
and  take  appropriate  action  before  disaster  occurs.  Monitoring  scour  depths  at  bridges 
serves  two  purposes.  First,  it  is  a  valuable  source  of  information  for  prototype  scale  data 
that  will  ultimately  lead  to  better  predictive  methods.  Second,  it  ascertains  whether  a 
scour  hole  is  in  danger  of  undermining  a  pile.  If  this  is  the  case,  scour  mitigation  methods 
must  be  implemented.  Having  identified  problem  areas,  the  engineer  must  devise  a 
schedule  for  executing  local  scour  count ermeasures.  The  engineer  would  be  greatly  aided 
in  this  task  by  knowing  how  much  time  he  had  to  implement  these  methods  for  each 
specific  site.  In  other  words,  a  knowledge  of  the  temporal  variation  of  local  scour  would 
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help  in  identifying  which  bridges  need  immediate  attention  and  which  simply  require 
continued  monitoring. 

While  the  prediction  of  the  equilibrium  scour  depth  has  received  considerable 
attention  by  numerous  researchers  world  wide,  the  temporal  variation  of  scour  has  been 
relatively  ignored.  The  importance  of  the  rate  of  local  scour  is  fourfold.  First,  any  study 
pertaining  to  the  mechanisms  of  local  scour  is  in  itself  beneficial  to  the  state  of  the  science. 
Second,  as  mentioned  previously,  the  rate  of  local  scour  can  help  in  monitoring  activities 
and  in  the  scheduling  of  scour  mitigation.  Third,  understanding  the  temporal  dependence 
of  local  scour  is  the  first  step  in  predicting  scour  under  conditions  of  unsteady  flows  such 
as  waves  or  tidal  currents.  Finally,  the  current  method  for  prediction  of  local  scour 
involves  using  the  design  flood  conditions  as  input  into  formulas  for  calculation  of  the 
equilibrium  scour  depth.  Quite  often,  the  time  taken  for  the  design  discharge  to  scour  out 
this  depth  is  larger  than  the  time  for  which  the  event  runs.  This  is  especially  true  for 
bridges  situated  in  coastal  zones.  Figure  1 . 1  shows  typical  examples  of  hydrographs  for 
both  floods  in  riverine  locations  and  storm  surges  associated  with  coastal  events.  The 
time  scales  of  both  events  are  vastly  different.  This  leads  to  the  reasonable  question  of 
whether  the  same  sets  of  equations  should  be  used  in  both  areas. 

1.2  Local  Scour 

When  a  pier  is  placed  in  an  alluvial  stream,  the  flow  pattern  around  the  obstruction 
is  modified  significantly.  Flow  is  accelerated  around  the  sides  of  the  pier.  Also,  an 
adverse  pressure  gradient  forms  upstream  of  the  pier.  The  upstream  boundary  layer 
undergoes  a  three  dimensional  separation  leading  to  the  formation  of  a  horseshoe  vortex. 
These  changes  in  the  hydrodynamics  significantly  increase  the  shear  stress  at  the  bed.  If 
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the  bed  is  erodible  and  the  shear  stresses  are  great  enough,  a  scour  hole  forms  around  the 
pier.  This  phenomena  is  known  as  local  scour.  Local  scour  can  occur  around  any 
obstruction  placed  in  a  current  over  an  erodible  bed  including  both  bridge  piles  and  bridge 
abutments. 


0  1  2  3  4  5 

Time  (days) 


Figure  1 . 1  Typical  hydrographs  of  a  flood  event  and  a  storm  surge  event. 

Local  scour  can  be  classified  as  either  clear  water  or  live  bed.  For  clear  water 
scour,  the  shear  stress  on  the  bed  away  from  the  pile  is  less  than  the  critical  shear  stress 
necessary  for  sediment  movement.  There  is  no  movement  of  sediment  on  the  bed  where 
the  influence  of  the  flow  modifications  by  the  pier  is  negligible.  In  this  condition,  the 
scour  hole  grows  rapidly  initially  and  then  slowly  approaches  an  equilibrium  value.  The 
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equilibrium  scour  condition  is  defined  as  the  depth  at  which  the  hydrodynamics  produced 
shear  stresses  are  no  longer  able  to  remove  sediment  from  the  hole. 

For  live  bed  scour,  sediment  is  in  motion  everywhere  on  the  bed.  The  shear  stress 
at  the  bed  away  fi-om  the  pile  is  greater  than  the  critical  value.  In  this  case,  scour  occurs 
at  a  much  greater  rate  than  the  clear  water  condition.  The  equilibrium  scour  depth  is 
attained  in  a  much  shorter  length  of  time.  The  equilibrium  condition  is  defined  as  the  point 
where  the  amount  of  sediment  entering  the  scour  hole  is  equal  to  the  amount  of  sediment 
being  removed.  Hence,  there  is  no  change  in  the  scour  depth. 

Chapter  2  contains  a  thorough  description  of  the  different  types  of  scour  as  well  as 
the  individual  mechanisms  of  scour  that  occur  throughout  the  process.  This  dissertation  is 
primarily  concerned  with  clear  water,  local  scour. 

1.3  Scope  of  Study 

The  objective  of  this  dissertation  is  to  provide  a  means  for  predicting  the  temporal 
variation  of  clear  water  local  scour  around  a  single  circular  cylinder.  In  the  following 
investigation,  laboratory  experiments  are  conducted  over  both  fixed  and  erodible  beds. 
The  commercially  available  software  FLOWS  D®  is  used  to  model  the  hydrodynamics 
associated  with  the  physical  experiments.  The  results  from  the  physical  experiments  are 
compared  with  those  from  the  numerical  simulations  to  validate  its  use.  Finally,  a  novel 
approach  to  predicting  the  rate  of  local  scour  is  developed  and  tested  using  results  from 
the  numerical  simulations. 

To  better  understand  the  nature  of  the  problem  of  local  scour.  Chapter  2  reviews 
the  current  state  of  knowledge  concerning  local  scour.  It  commences  with  a  brief 
summary  of  the  definition  of  local  scour  and  a  review  of  the  mechanisms  involved.  Next, 
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a  comprehensive  survey  of  the  literature  concerning  prediction  of  the  temporal  evolution 
of  local  scour  is  presented.  Finally,  the  research  involving  the  hydrodynamics  in  the  scour 
hole,  including  both  computational  and  experimental  investigations,  is  discussed. 

Chapter  3  reports  the  physical  experiments  conducted  as  part  of  this  investigation. 
Seven  physical  experiments  are  conducted  including  three  erodible  bed  experiments  and 
four  fixed  bed  tests.  From  the  experiments,  both  scour  and  hydrodynamic  data  are 
collected.  The  fixed  bed  tests  are  conducted  using  the  initial  flat  bed  condition  and  the 
scour  holes  formed  after  1,  6,  and  24  hours.  Dye  injection  studies  are  also  conducted  to 
gain  better  insight  into  the  hydrodynamics.  Plots  are  generated  showing  the  velocity 
vectors  measured  using  an  Acoustic  Doppler  Velocimeter  (ADV). 

Chapter  4  is  devoted  to  the  numerical  simulation  performed  using  the 
commercially  available  software  FL0W3D®.  FL0W3D®  is  a  three  dimensional  finite 
difference  model  that  solves  the  Navier  Stokes  equations  using  a  variety  of  turbulence 
closure  techniques.  Simulations  of  the  physical  experiments  in  Chapter  3  are  performed. 
An  additional  set  of  simulations  is  run  using  the  average  values  from  the  all  the  physical 
experiments  at  numerous  intermediate  scour  depths  between  the  initial  flat  bed  condition 
and  the  scour  hole  after  24  hours  of  scouring. 

Chapter  5  begins  with  a  comparison  of  the  results  from  the  physical  experiments 
and  the  numerical  simulations.  The  numerical  model's  performance  is  judged  on  its  ability 
to  predict  both  the  magnitude  and  direction  of  the  measured  velocity.  In  the  latter  part  of 
Chapter  5,  a  novel  approach  to  the  prediction  of  local  scour  is  developed.  This  method  is 
tested  using  the  second  set  of  simulations  in  Chapter  4  and  the  scour  data  from  Chapter  3 . 
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Yet  another  set  of  simulations  is  run,  modeling  a  completely  different  scour  situation  in 
order  to  validate  the  technique. 

Finally,  Chapter  6  presents  a  summary  of  the  work  performed  as  part  of  this 
dissertation.  Also  contained  in  this  chapter  are  recommendations  for  future  study. 


CHAPTER  2 
BACKGROUND  AND  LITERATURE  REVIEW 


In  this  chapter  a  review  is  presented  of  the  available  literature  relevant  to  this 
research.  The  first  section  includes  a  brief  survey  of  material  concerning  the  definition  of 
local  scour  around  bridge  piers  and  an  analysis  of  the  basic  mechanisms  involved.  The 
following  section  covers  the  available  papers  pertaining  to  prediction  of  the  temporal 
evolution  of  local  scour.  To  conclude,  research  involving  measurement  of  the  flow  field  in 
the  scour  hole,  including  computational  and  experimental  investigations,  is  discussed. 

2.1  Local  Scour 

Bed  elevation  changes  around  bridge  piers  caused  by  the  flowing  of  water  are  a 
major  reason  for  bridge  foundation  failure.  Bridge  piers  are  supported  by  fi-iction  between 
supporting  piles  and  the  surrounding  bed  material.  If  the  change  in  bed  elevation  is 
sufficient  to  uncover  these  piles,  the  total  friction  supporting  them  is  reduced  and  hence 
the  piles  may  settle  or,  in  some  cases,  the  bridge  may  fail  catastrophically.  Due  to  a  lack 
of  understanding  of  local  scour,  bridge  designers  have  responded  by  using  conservative 
predictive  equations.  In  order  to  reduce  the  costs  associated  with  the  construction  of  the 
foundations  while  maintaining  sufficient  safety,  more  information  must  be  gathered  about 
the  status  of  the  bed  elevation  under  design  event  conditions. 

Bed  elevation  change  can  be  attributed  to  a  number  of  different  factors  and  each 
must  be  taken  into  account  for  a  proper  design  .  These  changes  are  classified  under  four 
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different  categories.  The  first  two  types  occur  regardless  of  the  presence  of  the  bridge  and 
therefore  cannot  be  strictly  referred  to  as  bridge  scour.  The  last  two  are  directly  due  to 
the  construction  of  the  bridge  over  the  waterway.  First,  there  may  be  a  change  in  the 
overall  elevation  of  the  entire  waterway.  This  may  be  either  temporary  or  progressive  and 
either  involve  an  accumulation  or  loss  of  sediment  over  the  entire  bed.  These  phenomena 
are  known  as  aggradation,  in  the  case  of  sediment  accumulation,  or  degradation,  in  the 
case  of  sediment  loss.  Secondly,  there  may  occur  a  shift  of  the  thalweg.  The  depth  across 
a  waterway's  cross  section  may  not  be  uniform  due  to  the  presence  of  secondary  channels 
within  the  main  channel.  These  secondary  channels  may  shift  progressively  or  temporarily 
thus  uncovering  pile  foundations  that  were  not  previously  exposed.  Thirdly,  the  presence 
of  the  bridge  piles  reduces  the  cross-sectional  area  of  the  waterway  at  the  bridge.  This 
increases  the  average  velocity  of  the  flow  as  well  as  the  bed  shear  stress  over  the  section 
and  thus  increases  sediment  transport.  The  result  is  a  lowering  of  bed  elevation  across  the 
section.  This  is  known  as  general  or  contraction  scour.  The  final  type  of  scour,  local 
scour,  is  attributed  to  the  increase  in  shear  stress  in  the  immediate  vicinity  of  the  pile  due 
to  the  fact  that  the  presence  of  the  pile  alters  the  flow  field.  In  most  situations,  local  scour 
is  one  of  the  more  important  factors  in  determining  the  bed  elevation  at  the  individual 
piles. 

Johnson  and  McCuen  (1991,  p.  1)  define  local  scour  as  "the  erosive  action  of  the 
water  in  excavating  and  carrying  away  material  from  the  channel  bed.  An  obstruction, 
such  as  a  bridge,  causes  an  interference  in  the  flow  of  the  stream,  which  changes  the  flow 
pattern  at  the  obstruction.  This  results  in  a  deepening  of  the  scour  hole  around  the  bridge 
pier  beyond  the  level  that  would  naturally  occur  from  degradation  and  general  scour." 
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Figure  2.1  Shield's  diagram. 


In  practice,  local  scour  is  classified  by  whether  or  not  sediment  is  entering  the  hole 
from  upstream.  To  illustrate  this  classification,  the  critical  friction  velocity,  u«c,  is  defined 
as  the  friction  velocity  at  the  bed  required  to  initiate  sediment  motion  in  the  absence  of  the 
pier.  This  quantity  is  determined  empirically  through  the  use  of  Shield's  diagram  (Figure 
2. 1).  The  classification  is  based  on  the  value  of  the  ratio  of  the  bed  friction  velocity  to  the 
.  u» 

critical  friction  velocity,  — ,  where  the  bed  friction  velocity  is  defined  as  u»  = 

where  Xb  is  the  shear  stress  at  the  bed  and  p  is  the  mass  density  of  the  fluid.  The  subscript 
c  indicates  the  friction  velocity  corresponding  to  initiation  of  sediment  movement. 
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When  is  greater  than  1,  the  scour  is  said  to  be  "live  bed  scour."  Here, 

sediment  leaving  the  scour  hole  is  being  replenished  from  upstream  (if  sediment  is 
available),  though  at  a  lesser  rate.  In  the  area  away  from  the  pile,  the  shear  stress  at  the 
bed  is  greater  than  the  critical  value  and  hence  sediment  transport  is  occurring.  Under  live 
bed  conditions,  bed  features,  for  example  dunes,  may  be  present.  These  features  can 
travel  in  and  out  of  the  scour  hole  and  thus  under  "equilibrium"  conditions  the  scour  depth 
varies  temporally.  As  a  result,  equilibrium  is  defined  as  the  point  at  which  the  amount 
sediment  leaving  and  entering  the  scour  hole  are  equal. 

u* 

"Clear  water  scour"  occurs  when  1  >  >  0.5,  where  the  lower  limit  of  0.5  is  an 

approximate  value.  In  this  case,  sediment  is  removed  from  the  scour  hole  without  being 
replaced  by  a  supply  from  upstream.  In  the  area  upstream  of  the  pile,  the  sediment  is  not 
in  motion  and  the  friction  velocity  at  the  bed  is  less  than  the  critical  value.  Under 
conditions  of  clear  water  scour,  sediment  is  progressively  removed  from  the  immediate 
area  of  the  pile  until  the  flow  field  is  no  longer  able  to  remove  sediment  from  the  hole.  At 
this  point,  the  equilibrium  scour  depth  has  been  reached. 

u. 

For  completeness,  a  final  condition  is  defined  where  0.5  >  — .  There  is  no 

sediment  motion  near  the  pile  or  away  from  the  pile.  The  initiation  of  clear  water  local 

u, 

scour,  occurrmg  at  an  approximate  value  of  —  =  0.5,  has  been  verified  experimentally 
by  a  number  of  researchers  including  Hancu  (1971)  and  Ettema  (1980). 
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In  order  to  better  understand  the  scouring  process,  it  is  first  advantageous  to 
examine  the  individual  mechanisms  at  work.  These  processes  have  been  described  in 
detail  in  a  number  of  papers  (e.g.  Shen  et  al.  (1966)).  To  begin,  take  the  situation  where  a 
single  vertical  cylindrical  pile  is  placed  in  a  wide  stream  with  a  bed  consisting  of 
cohesionless  sediment.  Viewing  the  cylinder  from  above,  the  presence  of  the  cylinder 
causes  the  flow  to  accelerate  around  the  sides  of  the  pile  reaching  a  maximum  value  at 
approximately  90  degrees  to  the  flow  direction.  In  the  wake  of  the  cylinder,  the  flow  field 
consists  of  vertical  (wake)  vortices.  The  shedding  frequency  and  intensity  of  these 
vortices  is  dependent  upon  the  Reynolds  number  based  on  the  cylinder  diameter 
bU 

( R  =  — )  where  b  is  the  diameter  of  the  pile,  U  is  the  depth  averaged  velocity  and  n  is  the 
kinematic  viscosity. 


Figure  2.2  Flow  patterns  associated  with  local  scour. 
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Viewing  the  cylinder  from  the  side,  the  approach  flow  has  a  velocity  profile  that 
varies  with  depth  and  is  typically  logarithmic  in  shape  (Figure  2.2).  The  variation  of 
velocity  in  the  bottom  boundary  layer  causes  a  stagnation  pressure  gradient  to  develop  on 
the  leading  edge  of  the  cylinder.  Flow  is  thus  directed  downward,  inducing  a  jet  that 
impacts  the  bed.  Many  researchers  including  Shen  et  al.  (1965)  and  Chiew  (1984)  have 
identified  this  downflow  at  the  front  of  the  pier  as  the  main  cause  of  local  scour.  Once  the 
flow  reaches  the  bottom  it  curls  up  forming  a  vortex  at  the  base  of  the  cylinder.  The  ends 
of  the  vortex  are  swept  downstream  by  the  surrounding  flow  bending  the  vortex  into  a 
horseshoe  shape  when  viewed  from  above.  This  flow  feature  is  called  a  horseshoe  vortex. 
A  similar  phenomena  is  found  on  the  leading  edge  of  the  cylinder  at  the  water  surface. 
This  vortex  is  called  a  surface  roller. 

Scour  occurs  when  these  hydrodynamic  effects  produce  bed  shear  stresses  that 
exceed  the  critical  value  for  the  local  sediment.  During  initiation  of  scour,  sediment  is 
removed  from  the  sides  of  the  pile  where  the  flow  has  been  accelerated.  Ettema  (1980) 
and  Hjorth  (1975)  have  documented  that  the  points  of  initiation  of  scour  occur  at  roughly 
±45  degrees  from  the  direction  of  the  undisturbed  approach  flow.  Melville  (1975), 
however,  measured  the  maximum  shear  stress  at  the  unscoured  bed  at  points  ±100  degrees 
from  the  approach  flow. 

From  here,  the  scour  hole  works  its  way  to  the  front  of  the  pile  where  the 
impinging  jet  and  horseshoe  vortex  have  been  established.  As  the  jet  loosens  the  sediment, 
the  horseshoe  vortex  removes  it  and  the  scour  hole  begins  to  deepen.  As  more  sediment  is 
removed,  the  hole  expands  both  in  depth  and  in  width.  It  resembles  an  inverted  right  cone 
whose  axis  coincides  with  that  of  the  pile.  The  angle  that  the  slope  of  the  hole  makes  with 
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the  vertical  is  approximately  equal  to  the  dynamic  angle  of  repose  of  the  sediment.  There 
are  differing  opinions  as  to  how  the  horseshoe  vortex  behaves.  Melville  (1975)  observed 
that  the  horseshoe  vortex  is  initially  small  and  relatively  weak.  Once  scour  is  initiated,  he 
found  that  the  downflow  increases  as  more  of  the  fluid  attains  a  vertical  component  of 
velocity.  Due  to  this,  the  vortex  grows  rapidly  in  size  and  strength.  As  the  hole  enlarges, 
the  horseshoe  vortex  grows  in  diameter  and  descends  into  the  hole.  With  further 
deepening  of  the  hole,  Melville's  measurements  indicate  that  the  circulation  associated 
with  the  horseshoe  vortex  continues  to  increase  as  equilibrium  is  approached,  due  to  the 
increase  in  the  cross-sectional  area,  but  at  a  decreasing  rate. 

A  contradictory  view  of  this  behavior  was  put  forth  by  Baker  (1978)  who  deduced 
that  the  initial  circulation  associated  with  the  original  unscoured  condition  is  equal  to  that 
associated  with  the  equilibrium  condition.  In  other  words,  the  circulation  remains 
constant  throughout  the  entire  scouring  process.  The  assumptions  Baker  made  underlying 
this  claim  are  that  the  input  vorticity  to  the  horseshoe  vortex  system  remains  constant  and 
the  vorticity  diffusion  out  of  the  system  to  the  cylinder  increases  as  scouring  occurs  due 
largely  to  the  fact  that  scouring  increases  the  surface  area  of  the  cylinder.  These 
assumptions  were  based  on  dye  visualization  experiments.  Melville's  (1975)  hot  wire 
anemometer  measurements  indirectly  refute  this  observation.  He  made  measurements  of 
the  downflow  in  front  of  the  pier  that  show  that  the  velocity  of  the  fluid  entering  the  hole 
increases  as  the  hole  deepens. (Figure  2.3)  This  indicates  that  more  fluid  is  entering  the 
hole,  thus  the  assumption  that  the  circulation  remains  constant  may  not  hold  true. 

When  the  vortex  scours  out  sand  from  the  bottom  of  the  hole,  an  "avalanche"  of 
sediment  can  occur  from  farther  up  the  side  of  the  hole  causing  the  hole  to  widen.  Once  a 
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sediment  particle  is  put  into  suspension  by  the  horseshoe  vortex,  it  is  swept  around  the 
sides  of  the  cylinder.  Here,  a  number  of  things  can  happen  to  it.  If  it  interacts  with  the 
wake  vortices  it  can  be  brought  higher  in  the  water  column  due  to  the  upward  flow 
associated  with  the  center  of  the  vortex  and  then  transported  downstream  where  it  is 
eventually  deposited.  If  it  is  swept  directly  behind  the  cylinder,  it  experiences  the 
relatively  calm  wake  region  between  the  shedding  vortices.  Here  it  will  settle  out  to  form 
the  characteristic  mound  behind  the  pile.  During  certain  phases  of  the  scour  process, 
sediment  from  the  wake  mound  has  been  observed  to  reenter  the  scour  hole  in  the  form  of 


Figure  2.3  Pier  Downflow  Velocities,  Vmax,  measured  at  the  nose  of  a 
cylindrical  pier  in  air  flow:  Uo=12.2  m/s.  The  theoretical  distribution  of  V^ax 

according  to  Shen  et  al.  (1965)  is  also  shown  (dashed  line).  Data  from 
Melville  (1975)  with  a  flow  velocity  of  Uo=0.25m/s  and  pile  diameter  of  D  = 
50.8  mm  are  superimposed,  d,  /  D  =  1.8  (■),  0.94  (•),  and  0.0  (o,  flat  bed 
condition)  where  ds  is  the  scour  depth.  (Graph  taken  from  Chiew  (1984)). 


bedload. 
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The  wake  vortex  system  is  highly  dependent  on  the  cylinder  based  Reynolds 
number,  R.  Two-dimensional  wake  vortex  behavior  has  been  well  documented  by 
numerous  researchers,  including  Schlichting  (1955).  At  a  Reynolds  number  of  60  or 
below,  vortices  remain  attached  to  the  cylinder.  Above  this  value,  the  vortices  shed  in  an 
alternating  manner  forming  a  flow  structure  known  as  the  Von  Karman  vortex  street.  The 

fb 

frequency  of  the  shedding  is  controlled  by  the  Strouhal  number,  St  =  — ,  where  f  is  the 

shedding  frequency,  b  the  cylinder  diameter  and  U  the  depth  averaged  velocity.  The 
Strouhal  number  was  found  experimentally  to  be  solely  a  function  of  the  cylinder  Reynolds 
number.  At  extremely  high  Reynolds  numbers,  R  >  3.0  x  10^  the  vortex  system  is  no 
longer  coherent  and  the  wake  becomes  fully  turbulent. 

Melville  (1975)  made  flow  measurements  in  the  wake  region  using  a  hot  film 
anemometer.  His  measurements  indicate  a  departure  from  the  two-dimensional  behavior. 
The  cylinder  Reynolds  number  for  his  experiments  was  1.3  x  lO"*,  which  gives  a  Strouhal 
number  of  St  =  0.2.  This  number  translates  to  a  shedding  frequency  of  f  =  0.985  Hz.  He 
observed,  however,  spectral  peaks  in  his  measurements  between  1.13  and  1.17  Hz.  This 
translates  to  a  14  to  19%  increase  in  frequency.  He  also  observed  that  the  shedding 
frequency  decreased  as  the  depth  of  the  scour  hole  increased.  This  result  was  expected 
because  as  the  scour  hole  deepens,  the  mean  flow  velocity  will  decrease  at  that  point. 
Since  the  Strouhal  number  is  approximately  constant  in  this  range  of  Reynolds  number, 
this  will  manifest  itself  as  a  decrease  in  frequency. 
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The  depth  of  local  scour  is  a  function  of  a  number  of  parameters,  many  of  which 
are  interrelated.  Assuming  the  situation,  described  above,  of  a  single  circular  cylinder  in 
an  erodible  bed,  the  factors  influencing  local  scour  can  be  summarized  as  follows: 
Parameters  describing  the  fluid: 

p        =        the  fluid  density 

V        =        the  fluid  viscosity 
Note  that  both  p  and  v  are  functions  of  the  fluid  temperature  and  salinity. 
Parameters  describing  the  sediment: 

d        s        a  characteristic  grain  size 

(for  example,  the  median  grain  size,  dso,  or  equivalent  grain  size,  de) 
a        =        the  standard  deviation  of  the  grain  size  distribution 
Ps       =        the  density  of  the  sediment 
u»c      =        the  critical  friction  velocity 
Ws       =        the  fall  velocity  of  the  sediment 
(|>        s        the  dynamic  angle  of  repose  of  the  sediment 
Parameters  describing  the  flow: 

yo       =        the  depth  of  the  approach  flow 

U       =        the  depth  averaged  approach  flow  velocity 

k        =        the  roughness  of  the  approach  flow 

So       =        the  energy  slope 

Sb       =        the  bed  slope 

Xb       =        the  bed  shear  stress 
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Parameters  describing  the  cylinder: 

b        s       the  cylinder  diameter 

kc       =        the  cylinder  surface  roughness 
A  final  addition  to  this  group  of  variables  is  time,  t. 

Many  papers  in  the  literature,  for  example  Breusers  et  al.(1977)  and  Melville  and 
Sutherland  (1989),  have  performed  dimensional  analysis  on  these  variables  in  an  attempt 
to  identify  the  most  important  dimensionless  groups  affecting  local  scour.  Concerning 
equilibrium  scour  depth,  Ontowirjo  (1994)  examined  a  number  of  these  dimensionless 
groups.  He  analyzed  as  much  reliable  data  as  he  could  gather  from  both  the  literature  as 
well  as  in  house  investigations  to  determine  which  parameters  exhibited  the  strongest 
influence  on  the  equilibrium  scour  depth.  He  confirmed  the  findings  of  a  number  of  other 

U  yo 

researchers  in  identifying  the  three  most  influential  dimensionless  groups  — ,  — ,  and 

c 

where  Uc  is  the  depth  averaged  flow  associated  with  the  initiation  of  sediment 

b 

motion  upstream  of  the  pile. 

The  variation  of  the  equilibrium  local  scour  depth  with  the  depth  averaged  velocity 

ratio,  — — ,  has  been  thoroughly  examined  by  a  number  of  researchers.  It  can  be  shown 

u, 

that  this  parameter  is  proportional  to  the  ratio  of  friction  velocities,       ,  mentioned 

previously.  It  is  generally  agreed  that  the  initiation  of  scour  occurs  at  a  depth  averaged 
velocity  ratio  of  approximately  0.45  to  0.5.  From  this  value  the  equilibrium  depth  of  scour 
increases  almost  linearly  to  a  maximum  at  the  transition  between  clear  water  and  live  bed 
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scour,  — —  =  1.0.  In  the  live  bed  regime,  the  equilibrium  depth  of  scour  decreases  initially 

with  velocity  ratio.  It  then  approaches  a  second  maximum  as  the  bed  planes  out;  see 
Figure  2.4.  The  height  of  the  second  maximum,  relative  to  the  first,  is  still  a  subject  of 
some  debate.  Melville  (1984)  published  that  the  height  of  the  maximum  at  the  transition 
to  live  bed  was  a  function  of  whether  or  not  the  sediment  was  ripple  forming.  In  the 
analysis,  he  claims  all  the  data  collapses  to  either  a  ripple  forming  curve  (djo  <  0.6  to  0.7 
mm)  or  a  non  ripple  forming  curve  (dso  >  0.6  to  0.7  mm).  Between  both  curves,  the 
second  maximum  is  approximately  equal.  However,  in  the  ripple  forming  curve,  the  first 
maximum  is  less  than  the  second  and  the  reverse  is  true  for  the  non  ripple  forming  case. 
In  contrast  to  this  claim,  Sheppard  (1996)  believes  that  the  height  of  the 

maximums  are  a  function  of  the  sediment  to  structure  diameter  ratio,  —— .  His  analysis 

b 

shows  that  while  the  live  bed  maximum  appears  to  be  independent  of  grain  size,  the 
maximum  at  the  transition  to  live  bed  is  highly  sensitive  to  this  parameter.  He 
hypothesizes  that  rather  than  collapsing  to  a  ripple  or  non  ripple  forming  curve,  the  data  is 
represented  by  a  family  of  curves  in  between  these  two  extremes  that  are  dependent  upon 
sediment  to  structure  diameter  ratio,  see  figure  2.4.  The  considerable  scatter  that  is  found 
in  Melville's  data  may  be  explained  by  this  hypothesis.  Note  that  due  to  the  maximum  that 

occurs  in  the  dependence  of  —r—  (discussed  later,  see  Figure  2.6),  the  curves  shown  in 

b 

Figure  2.4  are  not  unique.  Rather,  the  curves  increase  from  the  lowest  to  the  highest  with 
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increasing       until  a  maximum  is  reached  and  then  decreases  back  down  to  the  lowest 
b 

curve  as  —  further  increased, 
b 

2.50  — , 
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Figure  2.4  Characteristic  dependence  of  non-dimensional 
equilibrium  scour  depth  on  depth  averaged  velocity  ratio. 


The  influence  of  the  aspect  ratio,  — ,  is  found  in  its  impact  on  the  surface  and 

b 

bottom  vortices  located  on  the  upstream  edge  of  the  pile.  As  discussed  earher,  the 
establishment  of  the  downflow  and  the  horseshoe  vortex  is  primarily  a  function  of  the 
bottom  boundary  layer  thickness.  Using  a  limiting  argument,  as  the  water  depth 
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approaches  infinity,  if  the  velocity  ratio  were  held  constant,  intuitively  there  should  come  a 
point  beyond  which  the  equilibrium  scour  depth  would  no  longer  be  a  function  of  water 
depth.  There  is  effectively  a  "practical"  boundary  layer  thickness  limit.  Laboratory  data 

show  that  if  —  and  —are  held  constant,  the  equilibrium  scour  depth  increases  rapidly 

V  y  0 

with  —  from  zero  and  asymptotically  approaches  a  constant  value  at  around       =  2.5  to 
b  D 

3  .0,  The  typical  behavior  of  scour  depth  with  aspect  ratio  is  shown  in  figure  2.5. 
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Figure  2.5  Characteristic  dependence  of  non-dimensional 
equilibrium  scour  depth  on  aspect  ratio. 
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The  importance  of  the  sediment  to  structure  diameter  ratio,  — ,  is  less  intuitive 

than  the  other  two  parameters.  To  begin,  the  range  of  laboratory  data  spans  several 
orders  of  magnitude.  For  this  reason,  it  is  convenient  to  work  with  the  logarithm  to  the 
base  10  of  this  parameter.  While  the  velocity  and  aspect  ratio  can  be  adequately  modeled 


the  laboratory,  the  sediment  to  structure  diameter  ratio  cannot.  As  a  result,  there  are  a 
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Figure  2.6  Characteristic  dependence  of  non-dimensional 
equilibrium  scour  depth  on  sediment  to  structure  diameter  ratio. 
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Figure  2.6  shows  the  typical  behavior  of  equiUbrium  scour  depth  with  log 


10 


\h  J 


where 


the  curve  in  the  prototype  scale  range  is  based  on  extrapolated  field  data.  It  is  interesting 
to  note  that  for  a  given  grain  size,  the  non-dimensional  equilibrium  scour  depth  first 
increases  then  decreases  with  increasing  structure  diameter. 

Why  this  parameter  correlates  so  well  with  the  laboratory  data  is  not  immediately 
obvious.  A  possible  explanation  put  forth  by  Sheppard  (1996)  is  that  this  parameter  is 
actually  the  ratio  of  two  different  Reynolds  numbers,  one  based  on  the  grain  diameter  and 
one  on  the  structure  diameter;  that  is, 


V 

where  the  velocity  very  near  to  the  bed,  u,  is  proportional  to  the  depth  averaged  velocity, 
U.  Both  of  these  Reynolds  numbers  characterize  the  flow  and  sediment  transport  in  the 
vicinity  of  the  pile. 


While  the  amount  of  literature  involving  equilibrium  scour  depth  is  plentifial,  the 
number  of  papers  that  pertain  to  the  time  rate  of  local  scour  are  relatively  scarce.  This 
section  will  first  summarize  a  number  of  observations  that  researchers  have  made  about 
the  temporal  variation  of  local  scour.  Next,  a  review  is  presented  of  the  models  available 
for  describing  this  variation. 


bU 


V 


(2.1) 


2.2  Temporal  Variation  of  Local  Scour 
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2.2.1  Observations  on  the  Time  Rate  of  Local  Scour 


The  equilibrium  scour  depth  is  approached  asymptotically  for  the  case  of  clear 
water  scour.  For  the  case  of  live  bed  scour,  the  equilibrium  scour  depth  is  approached 
more  rapidly  and,  once  reached,  the  depth  of  scour  oscillates  as  bed  forms  propagate  past 
the  pier.  The  height  of  these  oscillations  is  a  function  of  the  height  of  the  bed  forms.  Shen 
et  al.(1969)  reported  that  the  amplitude  of  the  fluctuations  is  less  than  Vi  the  height  of  the 
bed  forms.  The  distinction  in  temporal  development  between  the  two  is  illustrated  in 
Figure  2.7. 
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Figure  2.7  Scour  Depth  for  a  given  pier  and  sediment  size  as 
a  function  of  time  for  live  bed  and  clear  water  cases. 


Some  of  the  first  researchers  to  address  the  subject  of  local  scour  were  Shen  et  al. 
(1965).  In  this  study,  21  experiments  were  performed  using  a  single  cylinder  diameter  and 
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sediment  size,  but  varying  the  water  depth  and  the  depth  averaged  flow  velocity  so  as  to 
cover  both  clear  water  and  live  bed  conditions.  From  these  tests,  they  concluded  that  the 


depth  as  well  as  the  time  to  reach  equilibrium.  They  supplied  a  set  of  curves  for  use  in 
calculating  the  scour  depth  for  any  time  as  a  function  of  F.  The  fact  that  only  one 
sediment  size  and  cylinder  diameter  were  used  casts  some  doubts  onto  these  results. 

A  year  later,  Shen  et  al.  (1966)  attempted  a  more  fundamental  analysis  of  the  rate 
of  local  scour.  They  performed  37  experiments  involving  both  whole  and  half  cylinders  in 
fine  sediments  (dso  =  0.24  mm,  0.46mm)  over  a  range  of  live  bed  velocities.  They  also 
incorporated  the  work  of  Chabert  and  Engeldinger  (1956)  into  their  analysis.  Defining  t75 
as  the  time  it  takes  to  reach  75%  of  the  equilibrium  scour  depth,  they  reported  that  given 
similar  experimental  conditions  tys  tends  to  be  greater  for  larger  sediment  sizes.  They  also 
found  that  hs  was  smaller  for  larger  values  of  U^.  Both  of  these  results  make  intuitive 
sense.  They  also  reported,  however,  that  they  could  find  no  dependence  of  t75  on  the 
diameter  of  the  cylinder.  They  do  concede  however  that  the  range  of  cylinder  diameters 
tested  may  have  been  too  small  to  capture  any  systematic  variation. 

In  his  dissertation,  Melville  (1975)  made  a  detailed  description  of  the  formation  of 
the  scour  hole  under  clear  water  conditions  including  the  changes  in  the  hydrodynamics 
and  the  sediment  transport  as  a  function  of  time.  To  begin,  in  the  unscoured  condition, 
the  flow  ahead  of  the  transverse  diameter  of  the  pile  approximates  that  of  the  solution  for 
the  potential  flow  about  a  cylinder.  In  the  lee  of  the  cylinder  he  noted  a  strong  upflow  and 
a  weaker  backflow.  He  measured  an  increase  in  velocity  around  the  circumference  of  the 
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was  the  most  important  variable  in  determining  the  scour 
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cylinder  reaching  a  maximum  at  approximately  ±100°  measured  from  the  upstream  edge 
(approach  flow  direction).  These,  he  states,  are  the  points  where  scouring  commences. 

Indentations  form  at  these  points  of  maximum  velocity  due  to  the  increase  in  shear 
stress  at  the  bed.  The  flow  near  these  indentations  is  similar  to  that  in  the  lee  of  the  crests 
of  ripples.  As  a  result,  the  two  small  holes  quickly  work  their  way  around  the  perimeter 
towards  the  front  of  the  cylinder  and  meet  at  the  centerline.  In  the  case  of  ripple 
formation,  the  eroded  material  is  redeposited,  however  in  this  case  the  material  is 
transported  away  from  the  hole  downstream  by  the  flow.  Now,  a  small,  shallow  scour 
hole  is  present  that  stretches  around  the  perimeter  and  is  concentric  with  the  cylinder. 

Melville's  measurements  indicate  that  the  initial  horseshoe  vortex  and  associated 
downflow  are  relatively  weak.  With  the  growth  of  the  scour  hole,  more  fluid  gains  a 
downward  component  thus  increasing  the  strength  of  the  downwards  jet  and  hence  also 
the  strength  of  the  vortex.  The  jet  acts  to  erode  the  bed  at  the  upstream  edge  of  the 
cylinder.  The  eroded  material  is  carried  around  the  perimeter  of  the  cylinder  under  the 
action  of  both  the  accelerated  flow  and  the  horseshoe  vortex.  The  hole  grows  rapidly 
forming  the  shape  of  an  inverted  cone  that  is  coaxial  with  the  cylinder  and  whose  angle  is 
equal  to  the  dynamic  angle  of  repose  of  the  sediment.  The  horseshoe  vortex  descends  into 
the  hole  and  expands  to  fill  the  hole.  This  initial  phase  of  scouring  is  characterized  by  the 
increase  in  the  erosive  action  of  the  downflow  and  the  decrease  in  that  of  the  accelerated 
flow  about  the  sides. 

Once  the  scour  hole  becomes  large  enough  to  completely  contain  the  horseshoe 
vortex,  a  peaked  hump  forms  at  the  bottom  of  the  hole  concentric  with  the  cylinder.  This 
hump  is  built  up  by  the  upslope  flow  in  the  hole.  This  flow  is  very  strong  near  the  cylinder 
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and  decreases  radially  outward  against  the  force  of  gravity.  Shear  stresses  on  the  bed  in 
the  hole  vary  in  a  similar  manner.  Erosion  occurs  below  the  hump  where  sediment  is 
forced  upstream  and  into  the  flow.  This  material  is  then  carried  by  the  flow  into  the  wake 
region  where  it  is  deposited.  The  shear  stresses  above  this  hump  are  insufficient  to  cause 
erosion.  However,  the  removal  of  sediment  below  the  hump  increases  the  slope  angle  of 
the  hole  causing  the  sediment  to  "avalanche"  into  the  erosion  zone  below  the  hump  thus 
restoring  the  angle  to  the  dynamic  angle  of  repose  of  the  sediment.  Once  below  the  hump, 
this  sediment  now  experiences  the  erosive  action  of  the  downflow.  This  is  the  process  by 
which  the  hole  grows  in  depth  and  in  size. 

As  the  hole  deepens,  the  circulation  associated  with  the  horseshoe  vortex 
continues  to  increase,  due  to  the  fact  that  the  cross-sectional  area  of  the  hole  keeps 
expanding.  However,  the  rate  of  increase  in  circulation  diminishes  with  scour  depth.  This 
rate  is  controlled  by  the  supply  of  fluid  by  the  downflow  which  in  turn  is  controlled  by  the 
upstream  conditions  of  the  flow.  The  strength  of  the  downflow  near  the  bottom  of  the 
hole  decreases  as  the  hole  deepens  and  hence  the  erosive  power  decreases.  Also,  if  an 
armor  coat  forms  at  the  bottom  of  the  hole,  this  coat  will  act  to  diminish  the  rate  of 
scouring. 

Eventually,  equilibrium  is  reached.  This  occurs  when  the  time  mean  of  the  shear 
stress  at  the  bed  becomes  incapable  of  removing  any  more  sediment  from  the  erosion  zone 
behind  the  hump  in  the  scour  hole.  Melville  (1975,  p.  193)  describes  the  equilibrium 
condition  as  the  "condition  at  which  the  depth  of  scour  ahead  of  the  cylinder  is  just 
sufficient  so  that  the  magnitude  of  the  vertically  downwards  flow  ahead  of  the  cylinder  can 
no  longer  dislodge  surface  grains  at  the  bed."  He  surmises  that  since  the  magnitude  of  the 
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downflow  is  primarily  a  function  of  the  cylinder  diameter  and  the  magnitude  of  the 
approach  flow  velocity,  the  equilibrium  scour  depth  should  also  show  such  dependence. 
In  other  words,  the  approach  flow  depth,  yo,  only  has  an  indirect  influence. 

Although  equilibrium  is  reached  in  the  bottom  of  the  scour  hole,  erosion  continues 
to  occur  at  the  mound  downstream  of  the  cylinder.  The  mound  grows  progressively 
flatter  and  is  stretched  downstream  by  the  flow  out  of  the  hole.  This  flow  has  an  upward 
component  which  aids  in  removing  sediment  from  the  region  of  the  cylinder. 

For  clear  water  conditions  Ettema  (1980)  found  that  he  could  clearly  identify  three 
independent  stages  of  local  scour  about  a  circular  cylinder.  The  first  stage  is  termed  the 
initial  phase.  This  phase  is  characterized  by  the  transition  from  the  planar  bed  (d,  =  0)  to 
the  principal  erosion  phase.  Here  scour  is  controlled  at  first  by  the  acceleration  of  flow  at 
the  flanks  of  the  cylinder.  The  strong  downflow  at  the  leading  edge  of  the  cylinder 
develops  and  the  acceleration  at  the  flanks  rapidly  weakens.  The  downflow  and  associated 
horseshoe  vortex  then  take  over  and  control  the  next  stage  of  scour  which  Ettema  names 
the  eroding  phase.  During  this  phase,  the  horseshoe  vortex  rapidly  grows  in  size  until  the 
entire  vortex  is  contained  within  the  scour  hole.  The  final  phase  of  scour  is  the  equilibrium 
phase.  This  represents  the  stage  where  the  final  equilibrium  depth  has  been  achieved  for 
given  pier  size,  sediment  and  flow  conditions. 

Ettema  notes  that  each  of  these  phases  can  be  described  by  the  relationship 


where  both  Ki  and  K2  are  fianctions  of  the  variables  a,  — — ,  — — ,  and  .  This 

b     b  u^ 

relationship  appears  in  the  form  of  straight  lines  on  a  semilog  plot  of  the  depth  of  scour 
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time  history  (see  figure  2.8).  Ettema  does  not  indicate  the  nature  of  the  relationship 
between  Ki,  K2  and  the  dimensionless  variables  nor  does  he  specify  the  points  of 
demarcation  between  the  different  stages.  In  fact,  he  notes  that  more  work  is  required  in 
this  area. 


Time  (sec)  Time  (sec) 

Figure  2.8  Illustration  of  Ettema's  three  phases  of  scour 
(conditions  for  the  test:  dso  =  0.28mm,  yo  =  39.5cm,  b  =  1 1 .4cm, 
and  U  =  26.7cm/s,  data  from  Gosselin  and  Sheppard  (1995)). 

Finally,  Chiew  and  Melville  (1996)  examine  the  time  required  to  reach  equilibrium 
for  clear  water  scour  conditions.  Their  research  was  motivated  by  the  fact  that  a  number 
of  investigations  in  the  literature  were  reporting  equilibrium  scour  depths  and  drawing 
conclusions  from  these  values  on  tests  that  clearly  did  not  last  long  enough  to  reach 
equilibrium  values.  In  some  cases  these  tests  were  as  short  as  four  hours.  In  an  attempt 
to  standardize  the  criteria  for  reaching  an  equilibrium  value,  Chiew  and  Melville  collected 
data  from  35  experiments  that  covered  a  wide  range  of  cylinder  diameters,  flow  depths, 
and  approach  flow  velocities.  Two  different  sediment  diameters  in  the  coarse  range  were 
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employed.  The  tests  were  allowed  to  run  until  the  time  to  equilibrium,  te,  was  reached. 
The  authors  defined  te  as  the  time  when  the  scour  hole  develops  to  a  depth  of  at  which  the 
rate  of  increase  in  depth  does  not  exceed  five  percent  of  the  pier  diameter  in  the 
succeeding  24  hour  period.  This  criteria  yielded  values  for  te  as  high  as  three  days  for 
some  cases.  The  data  indicates  that  for  a  given  approach  flow  depth  and  velocity  ratio, 
the  time  to  equilibrium  increases  with  increasing  pier  diameter.  This  result  is  expected 
since  the  size  of  the  scour  hole  is  related  to  the  pier  width  and  a  larger  hole  requires  a 
longer  scour  time.  The  data  also  shows  that  te  increases  with  increasing  velocity  ratio, 
holding  the  other  variables  constant.  The  reasons  for  this  result  are  more  difficult  to 
explain.  A  greater  velocity  ratio  implies  a  greater  capacity  of  the  flow  to  remove  sediment 
and  hence  the  shorter  time  it  would  take  to  reach  the  equilibrium  condition.  However, 
greater  velocity  ratios  are  associated  with  greater  depths  of  scour  and  therefore  the  longer 
it  would  take  to  reach  equilibrium.  The  authors  suggest  that  the  latter  effect  is  dominant 
and  the  former  is  subordinate. 

From  their  data,  Chiew  and  Melville  derive  the  empirical  relation 

t,  =3xlO^-^e^^^"'"'^  (2.3) 
which  they  claim  is  applicable  for  values  of  0.5  <  (U  /  Uc)  <  1 .0.  The  authors  provide  as 
an  example  at  the  prototype  scale,  if  the  pier  diameter  is  2.0  meters  and  the  velocity  ratio 
is  1 .0  with  a  depth  averaged  velocity  of  2.0  meters  per  second,  it  would  take  34.7  days  to 
reach  the  equilibrium  depth.  This  kind  of  investigation  is  certainly  necessary  to  increase 
the  state  of  knowledge  concerning  scour  related  research.  The  authors,  however,  admit 
that  the  results  are  preliminary  and  only  based  on  a  limited  number  of  laboratory  scale 


31 

data.  It  is  also  a  concern  that  their  equations  only  indirectly  involve  the  influence  of  flow 

depth  and  sediment  size. 

2.2.2  Models  of  the  Time  Rate  of  Local  Scour 

Since  the  problem  of  determining  the  equilibrium  depth  of  local  scour  is  so 
difficult,  few  researchers  have  turned  their  attention  to  modeling  the  rate  at  which  scour 
occurs.  One  of  the  flrst  investigations  of  this  sort  was  performed  by  Carstens  (1966). 
One  of  the  objectives  of  his  study  was  to  formulate  similarity  criteria  for  the  temporal 
development  of  several  forms  of  clear  water  local  scour.  His  assumptions  include  the  fact 
that  in  areas  of  active  scour,  the  boundary  layer  thickness  is  negligible.  He  also  assumes 
that  the  velocity  distribution  in  the  vicinity  of  the  obstruction  is  a  function  of  the  geometry 
of  the  obstruction  and  the  scour  hole.  For  a  circular  cylinder  he  assumes  that  the  scour 
hole  at  any  time  has  the  shape  of  an  inverted  right  cone  concentric  with  the  cylinder  axis 
with  an  angle  equal  to  the  angle  of  repose  of  the  sediment,  ^.  Given  these  assumptions 
and  ignoring  the  effects  of  voids  ratio,  the  rate  of  sediment  transport  out  of  the  hole,  q,,  is 
equal  to  the  rate  of  change  of  the  volume  of  the  hole: 

dV 

q.  =  ^  (2-4) 

where  V  is  the  volume  of  the  scour  hole.  Putting  this  expression  in  terms  of  the  depth  of 
scour  yields: 
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Based  on  the  experiments  by  LeFeuvre  (1965),  Carstens  hypothesizes  that  the  non- 
dimensional  sediment  transport  function  has  the  form: 


UBd 


f(^(Ns  -  N^J,        obstruction  geometry,  particle  geometryj  (2.6) 


where  B  is  the  width  of  the  scour  hole,  L  is  a  characteristic  length  of  the  obstruction,  d  is  a 

characteristic  grain  size,  N,  is  the  sediment  number  defined  as  y  n — 7—7 ,  and  Nsc  is  the 

/V(s-l)gd 

sediment  number  corresponding  to  the  lowest  value  for  which  scour  will  occur.  Using  the 
data  of  Chabert  and  Engeldinger  (1956)  to  calibrate  his  model,  he  deduced  the  sediment 
transport  function  to  be 
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Combining  equations  (2.5)  and  (2.7)  and  solving  for  ds  gives  the  following: 
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Equation  (2.8)  is  an  equation  for  time,  t,  as  a  function  of  scour  depth,  ds.  This  method 
may  be  limited  by  the  fact  that  it  is  based  on  a  small  range  of  laboratory  scale  parameters 
and  does  not  take  into  account  the  effect  of  dso  /  b. 

Another  method  was  outlined  by  Shen  et  al.  (1966).  Having  realized  the 
deficiencies  of  their  previous  method  for  calculating  scour  depth  as  a  function  of  time, 
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mentioned  in  the  previous  section,  they  compiled  their  data  and  the  data  of  Chabert  and 
Engeldinger  (1956)  and  tried  to  fit  the  equation 

d3=K,  InCK^t)  (2.9) 
to  the  eroding  phase  of  scour  as  identified  by  Ettema  (1980).  They  tried  a  number  of 

combinations  of  variables  to  find  the  functional  relationship  of  Ki  and  K2  to  the  flow  and 

sediment  conditions.  The  resuh  is  the  following  equation: 


where 
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m  -  0.026e-^'^^yo  (2.11) 
and 


^Ut^ 


vyo^ 


vyo^ 


(2.12) 

Cunha  (1975)  notes  that  this  equation  seems  somewhat  overly  complicated.  He  then 
states  that  since  it  was  based  on  a  narrow  range  of  flow  and  sediment  conditions  it  is 
probably  not  very  well  suited  to  practical  application. 

Nakagawa  and  Suzuki  (1975)  buih  upon  the  work  of  Shen  et  al.  (1966,  1969)  to 
formulate  a  model  that  assumes  that  the  horseshoe  vortex  is  the  main  scouring  mechanism. 
In  an  analysis  that  is  similar  to  Shen's,  they  assume  that  the  scale  and  strength  of  the 
vortex  remain  constant  throughout  the  development  of  the  scour  hole.  Using  the 
definitional  sketch  in  Figure  2.9,  they  use  the  results  presented  by  Shen  et  al.  (1969)  for 
the  circulation  of  the  approach  flow,  Fi: 
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where  Uo  is  the  surface  velocity  of  the  approach  flow. 


(2.  13) 


Figure  2.9  Definitional  sketch  for  the  analysis  of  Nakagawa  and  Suzuki  (1975). 

The  circulation  in  the  triangle  HOL  is  given  by 

=  £aQade  =  -27CQa'  =  -27tav^,  (2.14) 
where  a  is  the  radius  of  rotation  of  the  vortex  and  Vv  is  the  tangential  velocity  at  the 
circumference  of  the  vortex  and  is  assumed  to  be  acting  at  the  particle  level.  Using  dye 
injection  experiments,  they  determined  a  to  be  0.177b  regardless  of  pier  Reynolds  number. 
Using  this  and  equating  fi  and  G2  they  calculated  an  expression  for  the  bottom  velocity  in 
terms  of  the  depth  averaged  approach  flow  velocity: 
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where  n  is  Manning's  n  and  k  \%  von  Karman's  constant.  Next,  using  the  definition  of  drag 
force  on  a  particle  they  presented  the  non-dimensional  shear  stress  at  the  bed  to  be 
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where  Co  is  the  drag  coefficient  of  the  particle.  The  sediment  transport  rate,  qb  can  be 
calculated  if,  given  that  pe  is  the  probability  that  a  particle  leaves  the  bed  per  second,  pd  is 
the  probability  that  a  particle  comes  to  rest  after  traveling  a  distance  X,  and  Ajd'^  is  the 
area  the  particle  occupies  at  the  bed  with  the  volume  of  that  particle  being  Aad^,  one 
integrates  the  equation  of  continuity  between  x  and  x  +  X,  holding  pd  constant,  where  x  is 
measured  in  the  direction  shown  in  Figure  2.9.  The  result  is  equation  (2. 17). 
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dp^  (if  motion  is  resticted  to  the  stagnation  plane) 
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They  then  used  the  sediment  transport  model  of  Paintal  (1971)  to  relate  0  to  pe.  They 
observed  experimentally  that  most  of  the  sediment  was  removed  from  behind  the  peaked 
hump  in  the  scour  hole.  The  radial  distance,  Kib  of  this  hump  to  the  cylinder  was  found  to 
be  0.25b,  and  this  distance  was  invariant  throughout  the  scouring  process.  It  was  next 
assumed  that  all  the  sediment  was  removed  fi^om  the  hole  via  this  area  and  that  this  area 
was  fed  sediment  by  shding  down  the  surrounding  slope.  The  area  of  the  sloped  region, 
Si,  can  be  represented  by  equation  (2.18). 

(l  +  2K)d,  d' 

S,  -  }  ^  +  7   (2.18) 

(l  +  K,)sin(l)    (l  +  KJbtan(t)sin(t) 


Therefore  the  rate  of  sediment  being  supplied  to  the  hole  is  SiCos(})(dds/dt).  So,  the 
continuity  equation  for  scour  reduced  to  a  one  dimensional  equation  in  the  x  direction  is 
equation  (2.19): 
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where  K  is  the  porosity  of  the  bed  divided  by  100.  Putting  in  the  derived  definitions  for 
each  of  these  variables  and  solving  using  the  initial  conditions  ds  =  0  at  time  t  =  0  gives  an 
expression  for  t  in  terms  of  ds  (equation  2.20): 

2  o 
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1-XJa,  b^^^     b  ^2(l  +  K,)tan(t)V  b  J  ^  3(1 +  K,)tan' (j) 
where  s  is  the  rate  which  the  particles  in  the  stagnation  plane  are  transported  outside  the 

hole.  In  their  comparisons  with  experimental  data,  it  appears  that  this  equation  does  a 

good  job  of  representing  scour  in  the  initial  phase  of  scour,  however,  once  the  rate  slows, 

the  equation  fails  to  track.  The  above  analysis  is  important  in  that  it  incorporates  the 

stochastic  nature  of  sediment  transport  into  the  scouring  process.  Its  failings  perhaps  lie  in 

the  number  of  assumptions  that  were  made  about  the  hydrodynamics. 

Hjorth  (1977)  also  uses  a  stochastic  model  to  describe  local  scour.  He  applies  the 

method  originally  proposed  by  Einstein  (1950)  to  illuminate  the  stochastic  nature  of  local 

scour.  This  model  also  assumes  that  the  shear  stress  at  the  bottom  of  the  scour  hole  is  a 

function  of  the  instantaneous  depth  of  the  hole  and  is  derived  from  vorticity  conservation 

considerations.  Hjorth  uses  the  reasoning  that  the  magnitude  of  the  vorticity  along  any 

vortex  tube  varies  inversely  with  the  cross-sectional  area  of  the  tube.  Thus,  if  a  vortex 

tube  has  a  cross-sectional  area  of  A  in  the  undisturbed  flow  and  an  area  of  A'  in  the 

leading  stagnation  plane  containing  the  horseshoe  vortex,  then  the  ratio  of  vorticities 

between  the  horseshoe  vortex  and  the  undisturbed  flow  is  A  /  A'.  He  then  reasons  that 

since  shear  stress  may  be  considered  proportional  to  vorticity,  the  same  ratio  applies. 

Using  the  fact  that  scour  commences  when  the  velocity  of  the  approach  flow  is  half  the 

critical  value  independently  of  pile  diameter,  Hjorth  deduces  that  A  =  4A'.  He  also 


(2.20) 
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assumes  that  as  the  hole  grows,  the  area  of  the  horseshoe  vortex  grows  from  A'  to  A'  + 
Ah,  where  Ah  is  the  cross-sectional  area  of  the  hole.  Therefore,  the  ratio  of  the  shear 
stress  in  the  bottom  of  the  scour  hole  to  that  of  the  undisturbed  flow  is  A'  /  (A'  +  Ah),  and 
this  ratio  decreases  as  the  hole  grows  in  size.  Assuming  that  the  hole  can  be  described  as 
an  inverted  cone  with  a  flat  section  of  width  0.25b,  the  cross-sectional  area  of  the  hole  is 


The  model  makes  use  of  a  stochastic  bed  load  model  that  assumes  that  the 
probability  of  particle  motion  for  both  the  approach  flow  and  for  the  bottom  of  the  scour 
hole  is  defined  as  the  probability  that  for  a  specific  particle  to  turbulent  eddy  contact  the 
erosive  strength  of  the  eddy  is  greater  than  the  resistive  force  of  the  particle.  The  model 
assumes  that  the  lift  on  the  particle  is  neghgible  and  the  distribution  of  entrainment  has  a 
Gaussian  distribution.  Both  of  these  assumptions  have  been  shown  by  Christensen  (1975) 
to  be  erroneous.  The  net  transport  rate  was  calculated  by  computer  at  each  time  step  as 
the  number  of  particles  ejected  out  of  the  scour  hole.  Hence  the  net  volumetric  change  of 
the  scour  hole  is  calculated  by  knowing  the  volume  of  the  particles  ejected  and  the 
porosity  of  the  bed.  When  compared  with  data  from  Chabert  and  Engeldinger  (1956),  the 
model  did  not  perform  well.  It  did  however  show  the  ability  to  reproduce  several  features 
of  local  scour  including  changes  of  time  scales  for  different  conditions  and  a  smooth 
approach  towards  an  equilibrium  value. 

As  of  late,  increasing  numbers  of  researchers  have  realized  the  importance  of  the 
temporal  variation  of  local  scour.  As  a  result  many  models  have  recently  appeared  in  the 
literature.  Yanmaz  and  Altinbilek  (1991)  use  an  approach  very  similar  in  many  regards  to 


(2.21) 
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Carstens  (1966).  Again,  assuming  the  shape  of  the  scour  hole  to  be  an  inverted  right 
cone,  the  rate  of  sediment  transport  out  of  the  scour  hole  defined  by  equation  (2.22): 
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At  this  point,  rather  than  use  Carstens'  sediment  transport  fiinction,  they  use  one 
published  later  by  LeFeuvre  et  al.  (1970).  Their  sediment  pickup  fiinction  has  the  form 
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where  Ue  is  the  volume  rate  of  bed  removal  per  unit  area  of  bed  per  unit  time,  u  is  the 
velocity  of  the  flow  at  the  particle  level,  and  Fi  is  defined  by  equation  (2.24), 
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where  a  is  the  bed  slope.  For  the  sediment  transport  out  of  the  scour  hole  the  following 
relationship  was  proposed: 
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where  f  is  a  coefficient  of  proportionality  which  is  a  Sanction  of  the  geometry  of  the  scour 
hole,  the  sediment  and  the  flow  properties.  To  determine  the  form  of  f,  the  authors 
performed  71  scour  tests  on  both  circular  and  square  piles  using  two  different  sediment 
sizes,  five  different  flow  conditions  and  3  different  pile  diameters  for  each  type  of  pile. 
Using  regression  analysis,  they  determined  the  following  form  for  f 
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Combining  equations  (2.22),  (2.25)  and  (2.26)  yields  the  differential  equation  (2.27)  which 
the  authors  solve  using  a  second  order  Runga-Kutta  method. 
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In  a  discussion  of  the  paper  published  later,  Yanmaz  and  Altinbilek  (1992)  admit  the 
short-comings  of  their  analysis.  The  model  seems  only  to  apply  in  the  range  of 
experimental  values  on  which  f  was  based.  This  is  primarily  due  to  the  sensitivity  of  the 
drag  coefficient.  To  illustrate,  in  their  tests,  the  authors  used  only  coarse  sediments  with  a 
typical  Cd  value  of  1.70  giving  Cd^^  =  4.42.  If  a  fine  sediment  (dso  =  0.2)  is  used,  Co  is 
approximately  6.43  and  hence  Cd^  **  =  183.2.  This  large  change  in  this  one  parameter  may 
cause  the  solution  to  "blow  up".  Also,  the  range  of  pier  diameters  is  also  quite  small: 
4.7cm  <  b  <  6.7cm.  Thus,  if  the  variables  used  are  outside  of  the  range  for  which  the 
model  was  calibrated,  unrealistic  resuhs  may  occur.  This  model  would  benefit  fi"om  the 
inclusion  of  a  wider  range  of  data. 

The  work  by  Johnson  and  McCuen  (1991),  was  motivated  by  the  desire  to  analyze 
the  scour  around  bridge  piers  caused  by  time  dependent  hydrologic  input.  Their  model 
consists  of  three  main  components  which  estimate  the  downflow  velocity,  the  termination 
velocity  (analogous  to  a  critical  velocity  for  the  unsecured  bed),  and  the  erosion  rate  and 
scour  depth.  They  assume  that  the  maximum  downflow  velocity  on  the  leading  edge  of 
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the  pier  is  a  function  of  the  cylinder  diameter,  depth  of  the  approach  flow,  and  the 
approach  velocity,  and  has  the  form: 

Va„ax=C,u(l-e-''"'»),  ■  (2.28) 
where  Ci  and  C2  are  constants.  The  form  of  equation  (2.28)  was  based  on  experimental 

observation.  This  maximum  velocity  occurs  at  some  unspecified  depth  along  the  leading 

edge.  The  downward  velocity  at  the  bottom  of  the  hole,  Vj,  is  found  using  equation 

(2.29) 

V,=KV,_,  (2.29) 
where  K  is  defined  by  equation  (2.30): 


K  =        .  .-C3d.-c.b  ■  (2.30) 
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C3  is  a  shape  parameter  and  C4  is  a  location  parameter.  The  termination  velocity  is  the 


velocity  where  the  downward  velocity  is  no  longer  capable  of  eroding  sediment  from  the 
bottom  of  the  hole.  In  other  words,  it  is  the  point  where  the  resistive  force  is  equal  to  the 
erosive  force.  The  erosive  force  in  this  case  is  the  shear  stress  at  the  bottom  of  the  hole. 
They  model  this  in  terms  of  Vd  in  equation  (2.3 1), 
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The  resistive  force  is  the  component  of  the  weight  along  the  slope  of  the  scour  hole,  or 
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where  Ps  is  the  submerged  mass  density.  Setting  equation  (2.32)  and  (2.3  l)(times  the 
cross-sectional  area  of  a  grain  (0.257rd^))  equal  to  each  other  and  solving  for  Va  (equal  to 
the  termination  velocity,  Vt,  at  the  critical  condition)  yields: 


2p  ed  sind) 

V,=J      p°       I      ,  (2.33) 
3pV2(l-sm(j)j 

where  the  subscript  u  denotes  the  fact  that  the  termination  velocity  has  not  been 
modified  to  account  for  turbulence  and  sediment  gradation.  Defining  the  sediment 
gradation,  G,  as  the  ratio  of  the  sediment  diameters  corresponding  to  the  84*  and  50* 
percentile,  d84  /  dso,  the  adjusted  termination  velocity  is  found: 


2pedsin(l) 

\=C,J  .  \y,  +C,Ge-"'",  (2.34) 

3pV2(l-sm(t)) 

where  Ce  accounts  for  the  inclusion  of  turbulence  and  C5  and  C7  are  fitted  coefficients. 
They  hypothesized  that  the  erosion  rate  was  proportional  to  the  difference  between  the 
downward  velocity  at  the  bed  and  the  termination  velocity.  They  proposed  the  following 
functional  form: 

-^  =  C,(V,-Vj.  (2.35) 
The  scour  depth  time  history  is  found  by  substituting  equations  (2.34)  and  (2.29)  into 
equation  (2.35),  and  then  choosing  an  appropriate  time  step  and  progressively  solving  for 
the  scour  depth,  or  vice- versa.  Constants  Ci  through  Cg  were  determined  using  a 
nonlinear  least  squares  method  to  fit  the  model  to  laboratory  data  from  Melville  (1975) 
and  Chiew  (1984).  In  using  this  model,  it  was  found  that  the  equations  were  extremely 
sensitive  to  time  step  size.  Also,  the  form  of  the  equation  for  erosion  rate  implies  that 
there  is  not  a  smooth  transition  to  the  equilibrium  value.  Once  the  downflow  velocity 
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drops  below  the  termination  velocity,  the  model  simply  stops  leading  to  a  discontinuity  in 
slope  once  equilibrium  is  reached. 

Sumer  et  al.  (1992)  published  an  ambitious  work  that  attempted  to  describe  the 
rate  of  scour  for  conditions  of  scour  under  both  steady  currents  and  periodic  flow  ranging 
in  time  scales  from  waves  to  tidal  flows.  For  the  case  of  steady  currents,  the  authors 
performed  1 8  experiments  all  under  live  bed  conditions.  They  propose  that  all  the  time 
history  curves  can  collapse  to  an  equation  of  the  form  given  in  equation  (2.36): 

d.-djl-e--),  (2.36) 
where  T  is  a  time  scale  that  is  a  function  of  the  flow  and  sediment  conditions.  The  form  of 

T  as  determined  from  the  experiments  is  given  by  equation  (2.37), 
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(g(s-l)d7''  2,000  b  lg(s-l)d. 
The  equilibrium  scour  depth  is  assumed  to  be  known  a  priori.  If  this  value  is  not  known, 

the  authors  recommend  a  value  of  1.0  -  1.5  times  b  based  on  the  results  of  Breusers  et  al. 

(1977).  Limitations  of  the  model  include  the  assumption  that  the  scour  time  history  can  be 

described  by  this  curve  in  all  situations  and  that  the  model  is  solely  based  on  data  in  the 

live  bed  regime. 

Kothyari  et  al.  (1992)  present  a  model  that  is  very  similar  in  methodology  to 
Hjorth  (1977).  They  begin  with  the  assumption  that  the  horseshoe  vortex  is  the  main 
scouring  mechanism.  For  the  initial  condition  of  the  unsecured  bed,  the  authors  use  an 
empirical  relation  to  find  the  diameter  of  the  horseshoe  vortex.  This  relationship  is  based 
on  experiments  performed  by  the  authors  as  well  as  those  performed  by  Qadar  (1980), 
Baker  (1979)  and  Muzzamil  et  al.  (1989).  This  relationship  has  the  form: 
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where  Dv  is  the  diameter  of  the  horseshoe  vortex.  They  also  assume  that  the  shape  of  the 
hole  is  an  inverted  cone  and  that  as  scour  progresses,  the  vortex  descends  into  the  scour 
hole  and  expands  to  fill  the  hole.  Hence,  the  cross-sectional  area  of  the  vortex  grows  from 
Ao  =  TcDv^  /  4  to  At  =  Ao  +As  =  (tcDv^  /  4)  +  (d^^/  2cot(l)),  where  At  is  the  total  area  of  the 
vortex  and  As  is  the  area  of  the  scour  hole.  Similar  to  Hjorth,  the  authors  assume  that  the 
shear  stress  under  the  vortex  is  a  function  of  its  cross-sectional  area.  Based  on  the  fact 
that  the  initiation  of  scour  occurs  at  approximately  one  half  the  critical  shear  velocity  of 
the  undisturbed  bed,  the  authors  propose  the  following  equation  for  shear  stress: 


Tp,  =4.0t 


(2.39) 


where  Tp,t  is  the  probabilistic  shear  stress  at  any  time  t  and  tu  is  the  shear  stress  on  the 
undisturbed  bed.  In  order  to  determine  the  time  scales  of  sediment  transport,  Kothyari  et 
al.  draw  upon  the  work  of  Paintal  (1971).  His  work  suggests  that  the  time  required  to 
remove  one  particle,  U,  is  given  by  equation  (2.40): 

C.d 

U=^^,  (2.40) 

where  u*,t  is  the  shear  velocity  based  on  Tp,t  and  po,t  is  the  average  probability  of  movement 
of  the  particle  at  time  t.  According  to  Paintal,  the  probability  is  given  by 
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The  constants  Ci  and  C2  were  determined  through  experiments  covering  a  wide  range  of 
sediment  and  cyhnder  diameters.  It  is  unclear  how  many  flow  conditions  were  considered. 
The  values  calculated  for  these  constants  were  0.57  and  0.05  respectively.  To  calculate 
the  scour  time  history,  the  authors  used  the  following  algorithm:  (1)  calculate  At  using  the 
depth  of  scour  and  equation  (2  .38),  (2)  calculate  the  shear  stress  at  the  bottom  of  the  hole 
using  equation  (2.39),  (3)  find  the  probability  of  movement  using  equation  (2.41),  (4)  find 
the  time  required  to  erode  one  particle  using  equation  (2.40)  and  add  this  to  the 
cumulative  time,  (5)  increment  the  depth  of  scour  by  one  sediment  diameter,  and  (6)  if  the 
shear  stress  calculated  in  step  (2)  is  greater  than  the  critical  shear  stress  return  to  step  (1). 
This  method  seems  quite  elegant  in  its  simplicity.  However,  it  is  unclear  how  the  time 
scale  can  be  based  on  the  removal  of  one  sediment  particle.  Also,  it  is  unlikely  that  the 
diameter  of  the  horseshoe  vortex  is  independent  of  the  velocity  of  the  flow  (see  equation 
(2.38)). 

The  models  described  above  cover  a  range  of  conditions  and  simply  give 
unreasonable  results  if  they  are  applied  outside  of  the  range  for  which  they  were  intended. 
As  such,  no  single  set  of  conditions  could  be  found  to  make  a  global  comparison  of  all  the 
models.  A  comparison  of  a  number  of  these  methods  was  made  by  Gosselin  and  Sheppard 
(1995).  They  chose  six  of  these  models  and  compared  them  to  experimental  data.  The 
results  are  given  in  Figure  2. 10.  The  inclusion  of  the  resuUs  fi"om  Sumer  et  al.  (1991)  is 
somewhat  unfair  since  their  model  was  intended  for  live  bed  conditions.  It  is  evident  from 
the  figure  that  the  prediction  of  the  rate  of  local  scour  is  still  an  investigation  worthy  of 
fiarther  research. 
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Figure  2.10  Comparison  of  predictions  of  rate  of  local  scour  with  experimental  data. 
(Test  conditions:  U  =  0.25m/s,  yo  =  0.11m,  dso  =  0.28mm,  b  =  0.051m) 


2.3  Experimental  and  Computational  Measurement  of  Local  Scour  Hydrodynamics 


Due  to  the  complexity  of  the  hydrodynamics  inside  the  scour  hole  as  well  as  the 


fact  that  the  bottom  boundary  changes  temporally,  it  is  very  difficult  to  make  flow 


measurements  near  to  the  pile  and  in  the  scour  hole.  As  a  result,  very  few  investigations 


have  been  performed  that  involve  measurements  of  the  flow  field.  One  of  the  earliest  j 

■■-1 

attempts  was  performed  by  Shen  et  al.  (1965).  His  analysis  of  temporal  variation  of  local  ; 


scour  included  measurement  of  the  flow  distribution  in  the  stagnation  plane  using  a  1/16* 
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inch  outside  diameter  Pitot  tube  connected  to  a  pressure  transducer.  In  order  to  make 
meaningful  flow  measurements,  the  erodible  bed  must  be  "fixed"  to  prevent  fiirther 
erosion.  The  investigators  used  two  methods:  (a)  a  water-glass-calcium  chloride  solution 
was  sprayed  on  the  bed  and  (b)  cement  powder  was  sprinkled  on  the  bed  and  allowed  to 
soak  in.  Neither  method  proved  to  be  exceptionally  effective,  however  it  was  possible  to 
collect  the  data  before  the  bed  degraded  significantly.  A  schematic  of  the  flow  patterns 
they  measured  is  contained  in  Figure  2.11. 
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Figure  2. 1 1  Schematic  of  flow  field  inside  the  scour  hole  from  Shen  et  al.  (1965). 


The  measurements  shown  involve  a  single  case  of  an  equilibrium  scour  condition. 
Their  analysis  yielded  a  number  of  results.  To  begin,  they  hypothesized  that  the  velocity 
profiles  at  the  bottom  of  the  hole  are  described  by  the  equation 
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u  z 

—  =  k  1-- 
U      ^  6; 


(2.42) 

where  u  is  the  vertical  velocity  as  a  function  of  the  distance  from  the  bottom,  z;  5  is  the 
thickness  of  the  boundary  layer  and  k  is  a  constant  of  proportionality.  The  measurements 
confirmed  this  result.  They  showed  that  the  separation  line  AB  in  Figure  2. 1 1  makes  an 
angle  of  7°  with  the  horizontal,  agreeing  with  the  concept  of  making  divergent  flow  not 
greater  than  7°  to  prevent  separation.  In  addition,  they  provided  vertical  velocity  profiles 
near  the  cylinder.  These  illustrate  an  increase  in  magnitude  of  the  velocity  from  y  =  '/2  b  to 
y  =  b  and  then  the  magnitudes  remain  constant  for  y  >  b  where  y  is  measured  from  the 
water  surface. 

The  following  year,  Shen  et  al.  (1966)  republished  some  of  their  findings  along 
with  a  number  of  new  observations  from  additional  experiments.  These  experiments  again 
included  measurements  of  the  flow  field  in  the  stagnation  plane  using  the  methods 
described  above.  They  chose  to  examine  three  different  scour  depths  over  a  range  of 
approach  flow  conditions  using  a  single  sediment  size  (djo  =  0.46mm)  and  two  different 
cylinder  diameters.  The  observations  made  included  the  fact  that  the  coming  flow  is  two- 
dimensional  and  that  it  follows  potential  flow  theory  at  least  in  the  stagnation  plane  near 
the  water  surface.  Potential  flow  theory  can  also  be  used  to  describe  the  pressure 
distribution  in  the  neighborhood  of  maximum  scour.  It  was  also  noted  that  the  boundary 
layer  thickness  in  the  scour  hole  was  constant  in  both  the  radial  and  the  tangential 
directions.  In  their  analysis  of  the  data,  Shen  et  al.  attempted  to  integrate  the  Navier- 
Stokes  equations  in  the  boundary  layer  at  the  bottom  of  the  hole  in  an  effort  to  verify  the 
boundary  layer  thickness  and  the  characteristic  velocity  as  a  function  of  the  flow 
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parameters  and  the  scour  depth.  The  semi-empirical  equation  that  resulted  from  this 
analysis  implies  a  dependence  of  the  boundary  layer  thickness  and  the  shear  stress  in  the 
scour  hole  on  the  Reynolds  number  to  the  -1/5*  power.  Unfortunately,  they  admit,  their 
analysis  is  flawed  in  that,  due  to  the  method  used  in  collecting  the  data,  measurement  of 
the  shear  stress  and  velocities  inside  the  boundary  layer  is  not  possible. 

As  part  of  his  dissertation,  Melville  (1975)  conducted  a  series  of  meticulous 
experiments  to  map  the  flow  field  at  three  different  stages  of  the  scour  process  for  a  single 
scour  event.  His  resuhs  were  also  included  in  a  paper  by  Melville  and  Raudkivi  (1977). 
The  purpose  of  this  investigation  was  to  provide  measurements  of  the  flow  field, 
turbulence  intensity,  and  boundary  shear  stress  in  the  scour  hole  as  well  as  the  interaction 
between  the  wake  vortices  and  the  horseshoe  vortex  in  order  to  better  understand  the 
individual  mechanisms  involved  in  the  scouring  process.  The  conditions  of  the 
experiments  were  as  follows:  yo  =  0. 1 5  m,  U  =  0.25  m/s,  dso  =  0.385  cm,  and  b  =  5.08  cm. 
The  stages  of  scour  chosen  were  at  t  =  0.0  hours  (no  scour  condition),  t  =  0.5  hours  and 
the  equilibrium  scour  condition.  The  bed  was  fixed  by  making  a  plaster  of  paris  cast  of  the 
scour  hole  and  then  using  this  to  cast  the  concrete  replica.  The  surface  of  the  replica  was 
coated  with  sand  used  from  the  experiment  to  maintain  the  bed  roughness.  Measurements 
of  the  flow  field  velocities  and  turbulence  intensities  were  made  using  a  unidirectional  hot 
film  anemometer.  The  anemometer  was  positioned  in  the  main  flow  direction  by  using  a 
short  tuft  of  cotton  glued  to  a  thin  rod.  In  addition,  Melville  used  the  hydrogen  bubble 
technique  to  trace  the  flow  patterns.  Measurements  were  made  in  four  planes  aligned 
radially  from  the  pier  at  the  angles  0°,  30°,  60°  and  90°  from  the  stagnation  plane  of  the 
cylinder.  In  addition,  measurements  were  made  in  the  wake  region  and  at  a  distance  of  2 
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mm  from  the  bed  throughout  the  scour  region.  The  shear  stress  was  estimated  by  using 
the  2  mm  velocity  in  the  equation 


Au 

T         =  K  

''0|2mm  ' 


(2.43) 


where  Au  is  the  velocity  measured  at  2  mm  from  the  bed  and  K  is  a  constant  of 
proportionality  which  was  determined,  through  experiments  in  the  approach  flow,  to  be 
equal  to  2.0.  Examples  of  these  measurements  are  found  in  Figures  2.12  and  2.13. 


Figure  2.12  Flow  direction  and  magnitude  in  the  stagnation 
plane  for  the  intermediate  scour  hole  from  Melville  (1975). 


A  number  of  conclusions  from  this  study  have  already  presented  in  this  chapter  in 
the  descriptions  of  the  scouring  process.  For  completeness,  they  will  be  repeated  below. 
Melville's  findings  can  be  divided  into  three  categories:  (a)  conclusions  based  on 
measurements  made  at  the  bed,  (b)  conclusions  resulting  from  measurements  in  the  four 
planes  around  the  cyhnder,  and  (c)  conclusions  from  measurements  in  the  wake  region. 


50 


Figure  2. 13  Flow  patterns  at  the  bed  for  three  stages  of  scour  (top:  flat 
bed,  middle:  intermediate,  and  bottom:  equilibrium)  from  Melville  (1975). 

Measurement  of  the  flow  at  the  bed  for  the  unscoured  condition  indicate  that 
scouring  is  initiated  by  high  local  shear  stresses  located  near  the  cylinder  at  approximately 
±100°  to  the  flow  direction  as  a  resuh  of  the  acceleration  of  the  fluid  around  the  cylinder. 
The  flow  pattern  at  the  bed  close  to  the  cylinder  is  significantly  altered  by  the  formation  of 
the  scour  hole.  The  shape  of  the  scour  hole  also  strongly  influences  contours  of  mean 
velocity,  RMS  velocity,  and  mean  shear  stresses  at  the  bed.  Mean  bed  shear  stresses  and 
turbulent  agitation  at  the  bed  tend  to  decrease  as  the  scour  hole  enlarges.  Also,  spectra  of 
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the  turbulent  velocity  fluctuations  near  the  bed  in  the  scour  hole  show  a  greater 
concentration  of  energy  in  the  1  to  10  Hz  range  than  is  found  in  the  approach  flow,  and 
also  lower  energy  at  higher  frequencies. 

Measurements  in  the  four  planes  surrounding  the  cylinder  indicate  that  a  strong 
vertically  downwards  flow  develops  at  the  front  of  a  cylinder  as  the  scour  hole  forms.  The 
velocity  at  the  bottom  of  the  hole  decreases  as  equilibrium  is  approached.  Melville 
hypothesizes  that  the  magnitude  of  the  downflow  velocity  is  directly  related  to  the  rate  of 
scour.  As  the  scour  hole  grows,  the  diameter  of  the  horseshoe  vortex  increases  and  the 
center  of  the  vortex  moves  away  from  the  cylinder.  The  circulation  associated  with  the 
vortex  increases  rapidly  initially  and  then  slows  and  approaches  a  constant  value  as 
equilibrium  is  reached. 

In  the  wake  region,  the  vortex  convection  speeds  decrease  with  depth.  This  is 
consistent  with  the  assumption  of  a  logarithmic  velocity  profile  bending  the  shed  vortices 
from  the  vertical.  In  addition,  the  frequency  of  shedding  showed  an  increase  from  the  two 
dimensional  case  by  approximately  14%  for  the  case  of  the  unsecured  condition.  The  data 
shows  a  decrease  in  frequency  as  the  scour  hole  enlarges.  This  result  is  expected  if  the 
Strouhal  number  is  to  be  maintained  at  a  constant  value  since  the  mean  flow  that  the 
cylinder  experiences  decreases  with  the  depth  of  the  hole.  Melville's  experiments  and 
subsequent  analysis  made  an  important  contribution  to  the  field  of  scour  research  in  that  it 
was  the  first  investigation  to  attempt  to  examine  the  mechanisms  of  scour  individually 
through  direct  measurement  and  record  the  behavior  of  these  mechanisms  throughout  the 
scour  process. 
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As  part  of  their  study  of  the  temporal  evolution  of  local  scour,  Kothyari  et  al. 
(1992  )  compiled  data  from  other  researchers  and  performed  their  own  experiments 
concerning  the  size  of  the  horseshoe  vortex  at  the  initial  stages  of  scour  (i.e.  at  the 
unsecured  condition).  They  hypothesized  that  the  thickness  of  the  boundary  layer 
upstream  of  the  pile  was  an  indicator  of  the  size  of  the  primary  vortex  of  the  system.  To 
measure  this  thickness,  they  injected  dye  into  the  flow  at  the  upstream  nose  of  the  cylinder 
using  a  3  mm  diameter  tube.  The  dye  followed  two  paths.  It  would  either  be  deflected 
sideways  and  be  carried  around  the  edge  of  the  pile  by  the  accelerated  flow,  or  it  would  be 
deflected  downwards  and  be  entrained  in  the  horseshoe  vortex.  The  data  gathered  using 
this  method  as  well  as  the  data  collected  from  Baker  (1979),  Qadar  (1980)  and  Muzzamil 
et  al.  (1989)  were  analyzed  to  produce  equation  (2  .44)  which  gives  the  diameter  of  the 
horseshoe  vortex,  Dv,  over  an  unsecured  bed  as  a  function  of  the  pier  diameter  and  the 
approach  flow  depth. 


—  =  0.28 


(2.  44) 


While  equation  (2.44)  does  show  the  intuitive  behavior  that  horseshoe  vortex  diameter 
increases  with  pier  diameter,  it  ignores  two  well  documented  features  of  the  vortex 
system.  First,  since  equilibrium  scour  depth  is  no  longer  a  function  of  aspect  ratio  beyond 
a  value  of  approximately  3.0,  it  would  then  follow  that  the  initial  vortex  diameter  should 
also  exhibit  this  asymptotic  behavior.  Second,  the  formation  of  the  horseshoe  vortex  is 
directly  related  to  the  vertical  variation  of  velocity  in  the  approach  flow.  This  velocity 
profile  is  a  function  of  the  approach  flow  roughness.  Hence,  it  seems  unlikely  that  the 
diameter  of  the  horseshoe  vortex  would  be  independent  of  sediment  size  or  at  least  bed 
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roughness.  In  addition,  the  authors  conclude  that  the  vortex  dimensions  are  independent 
of  approach  flow  velocity. 

Dey  et  al.  (1995)  published  an  ambitious  investigation  in  which  they  attempted  to 
describe  the  flow  field  around  a  circular  cylinder  in  an  equilibrium  scour  hole.  Analytical 
solutions  were  produced  for  the  flow  in  the  regions  in  the  scour  hole,  adjacent  to  the  pier 
above  the  flat  bed  and  in  the  wake  region.  The  authors  took  care  to  satisfy  the  continuity 
equations  and  matched  the  equations  using  empirical  constants  derived  from  experimental 
measurements  of  flow  velocities.  Their  experiments  employed  two  sand  diameters,  three 
pier  diameters,  three  approach  flow  depths,  and  six  approach  flow  velocities.  All 
experiments  were  performed  in  the  clear  water  regime.  The  length  of  the  scour  stage  of 
the  experiments  was  12  hours  when  the  authors  determined  a  quasi-equilibrium  state  had 
been  reached.  Then,  the  flume  was  drained  and  a  synthetic  resin  was  sprayed  on  the  bed 
to  stabilize  it.  Flow  measurements  were  made  using  a  five  hole  Pitot  sphere  with  an  outer 
diameter  of  8  mm.  In  areas  of  low  vertical  velocity,  a  three  tube  yaw  probe  with  an  outer 
diameter  of  3  mm  was  employed.  Measurements  were  made  in  the  planes  measuring  0°, 
15°,  30°,  45°,  60°,  and  75°  from  the  stagnation  plane  as  well  as  in  a  plane  perpendicular  to 
the  flow  (x  =  constant)  in  the  wake  region. 

Despite  the  number  of  simplifications  involved  in  producing  the  analytical 
solutions,  as  well  as  the  somewhat  dubious  method  of  data  collection  for  the  experiments 
used  to  produce  the  empirical  constants,  the  equations  for  the  flow  show  surprisingly  good 
agreement  with  measurements.  The  authors  provide  a  comparison  between  their  solutions 
and  the  measurements  performed  by  Melville  (1975)  reproduced  in  figure  2. 14.  The 
equations  themselves  are  somewhat  unwieldy  and  one  is  referred  to  the  paper  for  their 
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exact  form.  The  authors  should  be  lauded  for  their  novel  attempts  to  bridge  the  gap 
between  experimental  and  analytical  work  for  this  highly  complex  flow  field. 
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Figure  2. 14  Velocity  contours  in  the  stagnation  plane  (0°)  in  the  equilibrium  scour  hole 
from  (a)  Dey  et  al.  (1995)  (  predicted  from  equations)  and  (b)  Melville  (1975)  (measured). 
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Given  the  recent  advances  in  the  field  of  computational  fluid  dynamics,  it  is  logical 
that  researchers  would  utilize  this  very  powerful  tool  to  investigate  the  problem  of  local 
scour.  Olsen  and  Melaaen  (1993)  is  such  an  investigation.  They  coupled  non-transient 
three  dimensional  Navier-Stokes  equations  with  equations  for  bed  load  and  suspended 
load  for  the  sediment  transport.  The  Reynolds  stress  term  was  modeled  using  the  k-s 
method,  and  wall  laws  for  a  rough  boundary  were  used.  They  chose  to  model  an 
experiment  using  a  0.75  m  diameter  cylinder  in  a  bed  of  plastic  particles  measuring  3  mm 
in  diameter  with  a  relative  density  of  1 .04.  The  approach  flow  and  depth  were  0.067  m/s 
and  0.33  m  respectively.  The  experiment  was  stopped  after  12  hours,  well  before 
equilibrium  was  reached. 

The  goal  computationally  was  to  solve  for  the  hydrodynamics  using  the  Navier- 
Stokes  equations  and  then  use  this  solution  to  drive  the  sediment  transport  equations. 
Upon  solving  these  equations,  the  new  scoured  bed  was  determined  and  reintroduced  into 
the  flow  equations  to  find  the  new  flow  field.  This  process  was  repeated  until  the 
maximum  scour  depth  was  the  same  as  that  for  the  experiment.  The  equations  employed 
for  the  sediment  transport  were  a  convection-diffLision  equation  for  the  suspended  load 
and  Van  Rijn's  (1987)  transport  fianction  for  bed  load.  Based  on  sediment  continuity  for 
elements  close  to  the  bed,  erosion  and  deposition  of  sediment  were  calculated.  The  bed 
elevations  were  adjusted  so  as  to  allow  for  a  maximum  bed  elevation  change  of  10%  of 
the  approach  flow  depth.  Ten  iterations  between  the  flow  and  sediment  calculations  were 
required  using  this  procedure.  The  final  configuration  of  the  scour  hole  shows  good 

I 

agreement  with  that  of  the  experiment  (see  figure  2.15).  ] 
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Figure  2  .15  Bed  elevation  of  the  scour  hole.  Contour  line 
interval  =  2  cm.  From  Olsen  and  Melaaen  (1993). 


A  transient  solution  to  the  Navier-Stokes  equations  to  solve  the  hydrodynamics 
around  a  cylinder  was  developed  by  Richardson  et  al.  (1995).  They  used  the  commercially 
available  software  package  FL0W3D®  to  solve  the  three  dimensional  equations  for  the 
experimental  conditions  used  by  Melville  (1975)  for  the  cases  of  the  initially  flat  bed  and 
the  equilibrium  scour  hole.  They  chose  as  their  turbulence  closure  the  Prandtl  mixing 
length  model  in  order  to  reduce  computational  time.  They  found  that  the  software 
showed  relatively  good  agreement  between  computational  and  experimental  results.  They 
note,  however,  that  the  accuracy  of  the  solution  is  limited  by  the  choice  of  parameters 
which  affect  the  development  of  the  boundary  layer  and  caution  that  this  may  be 
problematic  when  coupled  with  a  sediment  transport  model.  They  go  on  to  recommend 
the  use  of  this  software  for  investigating  the  gross  flow  features  of  many  hydraulic 
situations. 


CHAPTER  3 
EXPERIMENTAL  APPROACH 

This  section  describes  in  detail  the  test  conditions  for  the  experiments  of  this 
study.  The  description  of  the  conditions,  equipment  and  experimental  methods  will  be 
comprehensive  so  that  the  results  may  be  properly  interpreted.  The  experiments  are 
divided  into  two  categories:  erodible  bed  experiments  and  fixed  bed  hydrodynamic  tests. 
Procedures  for  both  will  be  described  in  full.  This  chapter  will  conclude  with  a 
presentation  of  the  results  of  the  experiments  and  a  discussion  of  the  findings. 

3.1  Equipment 

The  equipment  used  in  these  experiments  fall  into  four  categories.  The  first  is  the 
flume  and  the  equipment  used  to  produce  and  measure  the  proper  flow.  The  second 
category  is  the  sediment  used  in  the  experiments.  The  third  category  is  the  model  used 
for  the  scour  experiments  and  the  method  for  fixing  the  bed  for  the  hydrodynamic 
experiments.  Finally,  the  equipment  used  for  the  measurement  of  the  flow  field  and  the 
dye  injection  studies  is  described. 
3.1.1  Flume 

All  of  the  experiments  were  performed  in  the  Hydraulics  Research  Flume  located 
in  the  Civil  Engineering  Department  of  the  University  of  Florida  (see  Figure  3.1).  The 
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main  section  of  this  recirculating  flume  has  a  flat  fiberglass  covered  bottom  and  is  30  m 
long,  2.46  m  wide  and  0.75  m  deep.  The  flow  is  produced  by  a  100  hp  pump  capable  of  a 
maximum  discharge  of  1 100  liters/sec.  The  volumetric  flow  rate  from  the  pump  is 
controlled  via  a  vertical  gate  (head  gate)  and  a  gate  controlled  30  inch  bypass.  The  flow 
then  enters  a  holding  area  which  is  taller  than  the  surrounding  flume  so  that  the 
volumetric  flow  rate  may  be  measured  via  a  V-notch  Thompson  weir.  The  weir  permits 
water  to  pass  to  the  flume's  main  section.  In  order  to  create  a  more  uniform,  consistent 
flow,  four  arrays  of  PVC  pipes  are  placed  in  hne  the  flow.  There  are  two  arrays  upstream 
and  downstream  of  the  weir.  The  arrays  consist  of  six  inch  long  pipes  that  are  2  inches  in 
diameter. 

Located  at  the  center  of  the  main  flume  is  a  windowed  test  section  which  is  0.33 
m  deeper  than  the  rest  of  the  main  flume.  This  section  starts  approximately  13  m  from 
the  weir  and  measures  6.10  m  long  by  2.46  m  wide.  This  section  is  filled  with  sediment 
such  that  the  level  of  the  sediment  equals  that  of  the  bottom  of  the  rest  of  the  main  flume. 
It  is  at  this  section  that  all  scour  experiments  and  hydrodynamic  measurements  are 
performed.  An  electric  motor  powered  carriage  rides  on  rails  mounted  on  the  top  of  the 
sides  of  the  flume  providing  a  stable  platform  on  which  observations  and  measurement 
may  be  performed  in  the  test  section.  At  the  end  of  the  main  flume  is  a  sluice  gate  (tail 
gate).  This  gate  permits  control  of  the  water  level  in  the  main  flume.  Behind  the  gate, 
flow  is  redirected  into  a  smaller,  deeper  return  channel  measuring  34  m  long,  1.22  m  wide 
and  1.00  m  deep.  The  return  channel  leads  back  to  a  reservoir  which  supplies  the  pump. 
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Figure  3.1  Schematic  of  Hydraulics  Research  Flume  (not  to  scale). 
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Flow  conditions  in  the  flume  are  controlled  by  the  head  gate  and  the  tail  gate,  with 
the  head  gate  controlling  the  volumetric  flow  rate  and  the  tail  gate  controlling  the  water 
depth.  Measurement  of  the  flow  depth  is  made  on  a  scale  mounted  to  the  windowed  test 
section.  Volumetric  flow  rate  is  measured  via  a  manometer  located  behind  the  V-hdtCh 
weir  in  the  holding  area.  The  value  of  the  discharge  over  the  weir  can  be  calculated  from 
the  head  behind  the  weir.  Using  BemouUi's  equation,  the  discharge,  Q,  is  determined  by 
the  equation 


8 


V2y 


Q  =  Cd-V2itan  -  H^'^  (3.1) 


where  0  is  the  spread  angle  of  the  weir,  H  is  the  head  behind  the  weir  and     is  an 
experimentally  determined  coefficient  which  is  a  function  of  Re,  We,  0  and  H.  Chadwick 
and  Morfett  (1993)  recommend  that  for  a  90°  weir  a  first  approximation  of  Cd  is  0.59. 
Using  this  value  of  C<j,  0  and  g  =  9.81  m/s^,  equation  (3.1)  becomes 


Q  =  1.394 


y  sec  y 


H" .  (3.  2) 


Prior  to  this  study,  the  discharge  over  the  weir  was  caUbrated  by  the  Civil 
Engineering  Department  at  the  University  of  Florida.  They  divided  the  surface  of  the 
flume  into  25  grid  locations.  At  each  grid  intersection,  the  flow  velocities  were  measured 
at  five  depths.  Analysis  of  this  data  resulted  in  the  plot  in  Figure  3.2  and  the  best  fit 
equation  given  in  equation  (3.3) 

Q  =  1594H"''*  (3.3) 
where  H  is  in  meters  and  Q  is  in  m^/sec.  This  relationship  is  used  here. 
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Figure  3.2  Calibration  of  V-notch  Thompson  weir. 


H  is  the  vertical  distance  from  the  bottom  of  the  weir  to  the  water  surface 
upstream  of  the  weir.  The  value  of  H  is  read  from  the  manometer  staildpipe  on  the  side 
of  the  flume  which  has  an  offset  of  5.62  cm  with  the  bottom  of  the  notch  in  the  weir. 
Therefore,  if  h*  is  the  value  read  from  the  manometer  in  centimeters,  the  value  of  H  to  be 
used  in  the  equation  is  (h*  -  5.62)/100.0. 
3.1.2  Sediment 

The  entire  test  section  was  filled  with  a  quartz  sediment  that  was  mined  in 
Jacksonville,  Florida.  The  properties  of  the  sediment  are  contained  in  Table  3.1  and  the 
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sediment  size  distribution  graphed  in  Figure  3.3.  The  grain  size  data  in  Table  3.1  and  tiie 
graph  in  Figure  3.3  were  produced  by  finding  the  best  fit  log-normal  distribution  curve 
using  the  mean  and  standard  deviation  as  the  free  parameters. 
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Figure  3.3  Sediment  size  distribution. 
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Table  3.1  Sediment  Properties 

Density,  p 

2650  kg/m' 

Median  Grain  Size,  dso 

0.169  mm 

Equivalent  Grain  Size,  de 

0.163  mm 

Geometric  Standard  Deviation,  Og 

1.30 

Uniformity  Coefficient,  Cu 

1.52 

Dynamic  Angle  of  Repose,  (j) 

31° 

63 

For  this  diameter  of  sediment,  ripples  will  form  on  the  bed  surface  given  a  stream  power 
per  unit  bed  area  between  approximately  0.09  and  1.1  Nm^s"'.  The  stream  power  per  unit 
bed  area  is  the  shear  stress  at  the  bed  times  the  depth  averaged  velocity. 
3.1.3  Scour  and  Hvdrodynamic  Models 

For  all  the  scour  and  hydrodynamic  tests,  a  vertical  circular  cylinder  with  uniform 
cross-section  was  used.  The  PVC  cylinder  was  1.0  m  long  and  16.83  cm  in  diameter.  For 
the  scouring  phase  of  the  experiments,  the  cyhnder  was  placed  in  the  test  section  such 
that  approximately  a  third  of  its  length  was  buried  in  the  sediment.  The  sediment  in  the 
test  section  was  then  leveled  to  a  height  equal  to  that  in  the  main  flume.  The  sediment 
was  also  smoothed  flat  over  the  entire  area  of  the  section.  The  reason  that  this  flat 
condition  was  chosen  as  a  starting  point,  rather  than  trying  to  mimic  the  bedforms  that 
might  be  found  in  a  prototype  channel,  was  to  have  a  consistent  starting  point  between 
experiments.  A  decal  with  quarter  centimeter  divisions  was  affixed  to  the  leading  edge  of 
the  cyhnder  (see  Figure  3.4(c)  and  (d)).  This  allows  for  visual  measurement  of  the  depth 
of  scour  at  the  leading  edge  during  the  scour  experiments. 

For  the  hydrodynamic  part  of  the  experiment,  the  bed  was  fixed  in  place  to 
prevent  changing  the  boundary  conditions  during  measurement  of  the  flow  field.  Once 
the  desired  scour  hole  had  been  produced,  the  flume  is  slowly  drained.  Then  a  single 
thickness  of  cheese  cloth  with  a  hole  cut  in  the  center  was  placed  over  the  cylinder  and 
onto  the  bed.  The  cloth  was  sprayed  with  a  mist  of  water  so  that  it  would  mold  to  the 
contours  of  the  scour  hole.  The  cloth  was  then  sprinkled  liberally  with  cement  dust.  The 
dust  was  then  sprayed  with  a  mist  of  water.  Finally,  sediment  of  the  same  type  used  in 
the  scour  experiment  was  sprinkled  over  the  top  of  the  cement  layer  to  provide  the  surface 
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with  the  same  roughness.  The  concrete  was  then  allowed  to  dry  before  the  flume  was 
filled  again  for  the  hydrodynamic  tests.  The  novel  aspect  of  fixing  the  bed  in  this  manner 
is  the  use  of  cheese  cloth.  The  cloth  strengthens  the  thin  concrete  layer.  It  was  found  that 
even  after  over  100  hours  of  hydrodynamic  testing,  there  was  virtually  no  degradiation  of 
the  bed.  Also,  the  labor  associated  with  this  technique  is  quite  low.  The  drawback  to  this 
technique  is  that  the  use  of  the  cheese  cloth  requires  much  more  cement  and  sediment  to 
cover  up  the  roughness  of  the  cloth.  As  a  result,  the  thickness  of  this  cheese  cloth  / 
sediment  /  cement  layer  is  greater  than  conventional  techniques.  For  these  experiments, 
the  thickness  of  the  layer  measured  approximately  3  mm.  Photographs  of  the  process  of 
fixing  the  bed  are  found  in  Figure  3.4. 
3.1.4  Flow  Measurement 

Point  measurements  of  the  flow  field  were  made  using  an  Acoustic  Doppler 
Velocimeter  (ADV)  manufactured  by  Sontek.  The  ADV  is  a  remote-sensing  three 
dimensional  velocity  sensor.  Utilizing  a  focal-point  acoustic  Doppler  system,  the  ADV  is 
inherently  free  of  signal  drift  and  therefore  does  not  require  constant  recalibration.  The 
system  used  in  these  experiments  is  the  ADVField.  This  system  consists  of  a  probe 
attached  to  a  7  mm  thick,  40  cm  long  stem  which  is  connected  to  a  signal  conditioning 
module  (see  Figures  3.5(a)  and  3.6).  The  module  is  connected  to  a  stand  alone  processor 
encased  in  a  "splash-proof  enclosure  via  a  10  m  cable.  The  processor  performs  the 
required  computations  for  real  time  estimations  of  the  3  components  of  velocities  at 
output  rates  up  to  25  Hz.  The  processor  also  acts  as  a  power  supply  for  the  probe. 
Communication  with  the  processor  was  made  through  the  serial  port  of  a  486  PC 
computer. 


Figure  3.4  "Fixing"  the  bed:  (a)  laying  down  the  cheese  cloth,  (b) 
sprinkling  with  cement,  (c)  and  (d)  before  and  after  photographs. 
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Figure  3.4  continued. 
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Software  provided  by  Sontek  allows  control  of  the  processor  and  also  enables 
storage  of  data  directly  to  the  hard  drive.  The  saved  file  is  a  binary  file  containing  all  the 
data  from  the  individual  experiment  and  information  can  be  extracted  from  this  file  via  a 
set  of  subroutines  provided  by  Sontek. 


Figure  3.5  The  Sontek  Acoustic  Doppler  Velocimeter  (ADV): 
(a)  Complete  system  and  (b)  close  up  of  measurement  probe. 


68 

Figure  3.5  (b)  shows  a  close  up  view  of  the  measurement  probe.  The  probe 
consists  of  a  center  transmit  transducer  and  three  receive  transducers  placed  around  the 
center  transducer  at  120°  angles.  The  operating  frequency  of  the  transmit  transducer  is  10 
MHz.  The  angle  between  each  receiver  and  the  transmitter-sampling  volume  center  line 
is  30°.  The  length  of  each  receiver  arm  is  27  mm.  The  distance  from  the  transmitter  to 
the  center  of  the  control  volume  is  cahbrated  by  the  manufacturer  for  each  probe.  For  the 
probe  used  in  these  experiments,  this  distance  is  53  mm. 

The  operation  of  the  system  begins  by  transmitting  short  acoustic  pulses  along  the 
transmit  beam.  As  the  pulses  move  through  the  water  they  are  continuously  being 
reflected  and  scattered  by  small  particles  found  in  the  water  (e.g.  micro-bubbles, 
suspended  sediment,  dust  particles,  seeding  material,  etc.).  The  three  receivers  are 
programmed  to  "listen"  to  the  signal  that  is  being  reflected  by  particles  found  inside  the 
sampling  volume  by  means  of  a  blanking  interval.  These  signals  are  Doppler  shifted  due 
to  the  relative  velocity  of  the  particles  in  the  flow  with  respect  to  the  probe.  It  is  assumed 
that  the  particles  in  the  flow  are  moving  with  the  same  velocity  as  the  fluid.  The  Doppler 
shift  detected  at  each  receiver  is  proportional  to  the  velocity  component  along  each 
bisector  between  the  send  and  receive  beanas. 

The  Doppler  shifts  thus  give  estimates  of  the  flow  velocity  along  three  different 
dkections  which  are  transformed  to  give  velocity  vectors  in  three  orthogonal  dkections. 
The  velocity  vector  is  centered  inside  the  sampling  volume.  This  volume  is  defined  by 
the  intersection  of  the  transmit  and  receive  beams.  It  has  the  shape  of  a  cylinder  with  a 
diameter  of  6  mm.  The  height  of  the  cyhnder  is  controlled  by  the  blanking  interval  and 
has  a  value  of  9  mm.  Thus  the  entire  sampling  volume  has  a  volume  of  approximately 
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0.25  cm^  and  the  value  of  the  velocity  vector  determined  by  the  processor  is  the  average 
value  over  the  whole  volume. 

Other  features  of  the  instrument  include  measurement  of  the  three  dimensional 
current  with  a  resolution  of  0.1  mm/s  over  the  range  of  0.0  to  2.5  m/s.  The  accuracy  is 
better  than  ±0.25%  or  ±0.25  cm/s  whichever  is  greater.  The  sampling  rate  is  conlxQllable 
up  to  a  value  of  25  Hz.  The  minimum  distance  from  the  sampling  volume  to  the  bottom 
boundary  is  5  mm.  The  operating  temperature  range  is  0°  to  40°C.  Much  of  the 
information  in  the  preceding  paragraphs  were  taken  from  the  articles  Kraus  et  al.  (1994) 
and  Lohrmann  et  al.  (1994).  Refer  to  these  articles  for  a  more  complete  description  of  the 
instrument  as  well  as  a  comparison  between  the  Sontek  ADV  and  Laser  Doppler 
Velocimeter  (LDV). 

One  of  the  naost  appealing  features  of  the  ADV  is  the  fact  that  there  is  a  minimum 
disturbance  of  the  flow  due  to  the  measurements  being  made  5  cm  away  from  the  probe 
tip.  To  minimize  the  flow  disturbance  even  further,  it  was  decided  that  the  probe  should 
be  mounted  in  such  a  way  as  to  allow  it  to  be  tilted  in  any  direction.  This  would  allow, 
for  example,  measurement  of  velocities  in  the  stagnation  plane  without  actually  placing 
the  instrument  anywhere  in  that  plane.  To  this  end  a  probe  support  allowing  this 
movement  was  designed  and  constructed  at  the  Coastal  and  Oceanographic  Engineering 
Laboratory  (see  Figure  3.6). 

The  mounting  bracket  allows  for  rotation  in  both  the  X-Z  and  Y-Z  planes  where 
positive  X  is  in  the  upstream  direction,  Y  is  across  the  flume,  and  positive  Z  is  vertically 
upwards.  The  rotation  bracket  which  holds  the  conditioning  module  in  place  runs  along  a 
vertical  track  (see  Figure  3.6).  The  bracket  is  moved  along  the  track  by  a  worm  gear  that 
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is  turned  by  a  stepper  motor.  The  stepper  motor  controller  allow  for  vertical  positioning 
to  millimeter  accuracy.  The  whole  assembly  is  mounted  to  a  traversing  mechanism  that 
consists  of  tracks  that  run  in  the  cross-tank  direction.  This  mechanism  also  uses  a  stepper 
motor  /  worm  gear  combination  for  movement  and  is  capable  of  millimeter  accuracy  in 
positioning.  The  traversing  mechanism  is  mounted  to  the  carriage  which  runs  in  the  long- 
tank  direction  and  is  positioned  via  a  scale  mounted  to  the  side  of  the  flume.  Positioning 
for  the  carriage  is  accurate  to  ±2  mm.  Thus  accurate  repositioning  of  the  probe  is 
possible  in  all  three  directions. 


Figure  3.6  Sontek  ADV  mounted  to  carriage:  (a)  in  vertical  and  (b)  in  tilted  position. 

In  order  to  calculate  the  rotation  of  the  probe  in  both  directions,  a  tilt  sensor  was 
tnountied  to  the  top  of  the  rotation  bracket.  It  is  the  black  box  mounted  at  the  top  of  the 
conditioning  module  in  Figure  3.6.  The  tilt  sensor  used  is  the  Model  900  Biaxial 
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Clinometer  manufactured  by  Applied  Geomechanics  Inc.  The  sensor  is  capable  of  +45 
degrees  of  rotation  and  has  a  resolution  of  0.02  degrees  of  an  arc.  The  output  of  the 
sensor  is  read  as  a  DC  voltage  between  ±2.0  volts.  The  cahbration  of  the  sensor  was 
performed  by  the  manufacturer.  To  convert  the  voltage  to  an  angle  of  rotation,  the 
manufacturer  provides  both  a  look-up  table  and  a  cahbration  curve  in  the  form  of  a  fifth 
order  polynomial.  The  curves  are  represented  graphically  in  Figure  3.7. 
3.1.5  Flow  Visualization 

The  flow  visualization  tests  were  performed  using  dye  injection  into  the  flow. 
The  dye  was  water  mixed  with  red  food  coloring.  It  was  injected  using  a  syringe  attached 
to  a  flexible  tube  with  an  inside  diameter  of  3  mm.  The  flexible  tube  was  wrapped 
around  a  metal  rod  shaped  into  an  "L"  so  that  it  could  be  held  in  position.  The  lighting 
for  the  visualization  was  provided  by  a  single  650  watt  lamp.  To  minimize  reflection  off 
the  water  surface,  the  lamp  was  beamed  through  the  bottom  of  a  clear  plastic  bucket. 
Photographs  of  the  dye  were  made  both  through  the  glass  window  and  through  the  water 
surface.  Examples  of  these  photographs  are  shown  in  section  3.2.3. 

3.2  Experimental  Methods  and  Conditions 
This  section  contains  a  comprehensive  review  of  the  experimental  part  of 
this  investigation.  Although  they  were  conducted  concurrently,  the  section  is  divided  into 
two  parts:  erodible  bed  experiments  and  fixed  bed  hydrodynamic  tests.  The  purpose  of 
the  erodible  bed  experiments  was  to  (1)  determine  the  time  history  of  depth  of  scour,  (2) 
find  the  equihbrium  scour  depth,  (3)  select  two  intermediate  conditions  for  measurement 
during  the  fixed  bed  hydrodynamic  tests  and  (4)  develop  the  scour  hole  to  be  fixed  in 
order  to  perform  the  flow  measurements.  In  total,  three  erodible  bed  experiments  were 
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performed.  Four  fixed  bed  experiments  were  performed  including  ones  with  the  same 
scour  depth  as  three  holes  discussed  above  plus  the  initial  flat  bed  condition.  The 
chronological  sequence  of  the  tests  are  as  follows:  (1)  flat  bed  hydrodynamic  test,  (2) 
scour  test  to  determine  equilibrium  scour  hole  (ds  =  13.8  cm),  (i)  equilibrium  hole 
hydrodynamic  test,  (4)  scour  test  stopped  after  one  hour  (ds  =  6.3  cm,  approximately  46% 
of  the  T  =  24  hour  depth),  (5)  one  hour  hole  hydrodynamic  test,  (6)  scour  test  stopped 
after  six  hours  (ds  =  10.5  cm,  approximately  76%  of  the  T  =  24  hour  depth),  and  (7)  six 
hour  hole  hydrodynamic  test. 
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Figure  3.7  Tilt  sensor  calibration  curves. 
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3.2.1  Erodible  Bed  Experiments 

The  erodible  bed  tests  involve  clear  water  scour  about  a  circular  cylinder.  All 
scour  experiments  were  performed  using  the  same  pile  diameter  and  sediment  size.  As 
mentioned  previously,  there  is  some  variability  in  the  flow  conditiOHS  both  between  and 
during  tests.  A  summary  of  the  average  flow  conditions  for  the  three  scour  tests  as  well 
as  the  other  relevant  parameters  is  presented  in  Table  3.2. 

Before  initiation  of  a  scour  test,  the  sediment  in  the  test  section  is  compacted  and 
leveled  to  the  same  elevation  as  the  rest  of  the  main  flume.  The  flume  is  then  slowly 
filled  with  water  to  the  appropriate  depth  and  allowed  to  sit  over  night.  When  begmning 
an  experiment,  the  sluice  gate  is  lowered  shut  and  the  pump  is  then  turned  on.  As  the 
discharge  quickly  ramps  up,  the  sluice  gate  is  progressively  raised  until  the  desired 
conditions  are  met.  This  method  allows  the  conditions  to  be  approached  from  the  lower 
velocity,  higher  depth  side  so  that  a  live  bed  condition  will  be  avoided.  The  approximate 
time  to  achieve  the  target  conditions  is  10  to  15  minutes.  As  the  flow  conditions  are 
being  reached,  the  depth  of  scour  at  the  front  of  the  pile  is  being  monitored.  The  defining 
criterion  for  the  start  of  the  experiment  is  the  point  at  which  scour  initiates  at  the  sides  of 
the  cyhnder.  Since  scour  occurs  rapidly  at  the  initial  stages,  measurements  of  the  depth 
of  scour  were  made  every  5  minutes  at  first.  This  interval  was  gradually  increased 
reaching  a  maximum  of  one  hour  at  the  end  of  the  equilibrium  scour  test.  The  water 
depth,  manometer  value  and  water  temperature  were  recorded  every  hour  for  the 
equilibrium  scour  test,  15  minutes  for  the  one  hour  test  and  half  hour  for  the  six  hour  test. 
At  the  end  of  a  test,  the  sluice  gate  was  again  lowered  thus  raising  the  water  depth  and 
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lowering  the  average  velocity.  Then  the  pump  was  turned  off  and  the  flume  was  allowed 
to  drain  slowly. 


Table  3.2  Summary  of  Experimental  Conditions  (Scour  Tests) 

Equilibrium 

1  Hour  Test 

6  Hour  Test 

Water  Depth,  yo  (m) 

0.348 

0.345 

0.349 

Discharge,  Q  (m  /s) 

0.210 

0.209 

0.211 

Depth  Avg.  Velocity,  U  (m/s) 

0.244 

0.245 

0.245 

Water  Temperature  (°C) 

24.7 

27.0 

27.7 

Final  Scour  Depth,  ds  (cm) 

13.8 

6.3 

10.5 

Cylinder  Diameter,  h  (cm) 

16.9 

Median  Grain  Size,  dso  (mm) 

0.169 

The  duration  of  the  test  to  determine  the  equilibrium  scour  depth  was  24  hours.  It 
was  decided  that  the  conditions  for  achieving  equilibrium  were  no  change  in  the  depth  of 
scour  after  four  hours.  The  scour  depth  at  the  leading  edge  after  24  hours  was  almost  14 
cm.  Although  it  is  unlikely  that  this  is  in  fact  the  true  equilihrium  depth  of  scour,  the  test 
was  stopped  at  this  point  because  the  time  and  cost  it  would  take  to  achieve  the  actual 
equilibrium  scour  hole  would  be  excessive.  In  fact,  according  to  the  equations  given  by 
Chiew  and  Melville  (1995),  for  the  conditions  of  this  experiment,  it  would  take  232  hours 
to  reach  equilibrium.  Using  their  curves,  after  24  hours  the  scour  depth  is  approximately 
80%  of  the  equilibrium  value.  The  time  history  of  the  depth  of  scour  is  shown  as  the 
sohd  line  in  Figure  3.8. 

One  of  the  purposes  of  the  equihbrium  test  was  to  identify  two  intermediate 
values  for  study  during  the  fixed  bed  tests.  By  plotting  the  depth  of  scour  time  history 
versus  the  logarithm  base  10  of  time  (Figure  3.9),  the  points  of  T  =  1  hour  and  T  =  6 
hours  were  chosen  such  that  the  former  would  occur  in  the  initial  rapidly  occurring  phase 
of  scour  and  the  latter  would  occur  after  the  rate  had  decreased  and  started  the  slow 
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approach  toward  an  equilibrium  value.  The  time  histories  for  both  of  these  tests  is  found 
in  Figure  3.8. 
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Figure  3.8  Depth  of  scour  time  history  for  the  three  erodible  bed  experiments. 
Once  the  tests  were  completed,  the  flume  was  drained  and  the  sediment  in  the  test 
section  was  allowed  to  dry.  The  hole  was  then  photographed  and  fixed  for  the 
hydrodynamic  tests.  Upon  completion  of  the  hydrodynamic  tests,  the  hole  was  mapped 
out  using  a  point  gauge  attached  to  the  traversing  mechanism.  Photographs  of  the  scour 
holes  are  found  in  Figure  3.10  and  contour  plots  of  the  scour  holes  are  found  in  Figure 
3.11. 
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Figure  3.9  Semilog  plot  of  depth  of  scour  time  history. 


3.2.2  Fixed  Bed  Hydrodynamic  Tests 

Four  fixed  bed  conditions  were  examined  for  this  part  of  the  investigation:  the 
three  scour  holes  discussed  previously  and  the  initial  flat  bed  condition.  Due  to  the 
number  of  point  measurements  of  the  flow  field,  the  hydrodynamic  testing  lasted  several 
weeks  for  each  case.  For  the  duration  of  the  testing,  the  flow  conditions  were  recorded  at 
half  hour  intervals-  The  average  flow  conditions  for  each  of  the  experiments  are  found  in 
Table  3.3.  Since  the  temperature  is  required  as  input  for  the  calculation  of  the  velocity 
vector  by  the  Sontek,  it  was  recorded  every  10  minutes. 


Figure  3.10  Photographs  of  the  scour  holes  after:  (a)  24  hours,  (b)  1  hour  and  (c)  6  hours. 
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Scour  Hole  After  24  Hours 


50.00 


Figure  3.1 1  Contour  plots  of  the  scour  holes  for  T  =  (a)  24  hours, 
(b)  1  hour,  (c)  6  hours.  Flow  direction  is  from  right  to  left. 
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Table  3.3  Summary  of  Experimental  Conditions  (Hydrodynamic  Tests) 

Flat  Bed  Test 

Equilibrium  Test 

1  Hour  Test 

6  Hour 

yo  (m) 

0.348 

0.348 

0.350 

0.348 

Q  (m'/s) 

0.211 

0.208 

0.210 

0.210 

U  (m/s) 

0.246 

0.242 

0.243 

0.246 

ds  (cm) 

0.0 

13J 

10.5 

6.3 

The  location  of  the  measurement  points  for  the  hydrodynamic  tests  varied.  For 
the  flat  bed  case,  the  first  case  studied,  a  much  denser  mesh  was  used.  The  purpose  for 
this  was  to  try  to  capture  the  highly  dynamic  flow  field  that  is  present  at  the  initial  stages 
of  scour.  In  addition,  measurements  were  made  in  two  planes  perpendicular  to  the  flow 
in  the  lee  of  the  cylinder  at  one  and  two  cyhnder  diameters  behind  the  center  of  the 
cylinder  to  try  to  capture  the  character  of  the  flow  field  in  the  wake  region.  Since  the 
flow  is  assumed  to  be  symmetric,  only  one  side  of  the  cylinder  was  studied.  This 
significantly  reduced  the  number  of  measurement  points.  A  graph  of  the  (X,Y) 
coordinates  of  the  measurements  for  this  test  is  presented  in  Figure  3.12(a). 

For  all  the  hydrodynamic  tests  the  following  coordinate  system  is  used:  positive  X 
points  into  the  flow,  positive  Z  is  upwards  and  positive  Y  is  such  that  the  right  hand 
coordinate  system  is  formed  (to  the  left  when  viewed  from  above  facing  upstream).  The 
location  of  the  origin  is  at  the  center  of  the  cylinder  at  the  level  of  the  unscoured  bed.  For 
the  remaining  three  tests,  measurements  were  made  in  planes  at  0°,  30°,  60°  and  90°  as 
measured  from  the  stagnation  plane.  A  graph  of  the  (X,Y)  coordinates  of  the 
measurement  positions  for  these  three  tests  is  presented  in  Figure  3.12(b). 
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Figure  3.12  Measurement  locations  for  (a)  flat  bed  case  and  (b)  remaining  cases. 
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The  vertical  location  of  the  points  was  identical  for  all  the  tests.  The  value  of  the 
Z  coordinates  were  as  follows:  24.0,  20.0, 15.0, 10.0, 7.5,  5.0, 4.0,  3.0,  2.0,  and  1.0  cm. 
In  addition  to  these  points,  measurements  below  the  unscoured  bed  level  (i.e.  in  the  scour 
hole)  were  made  at  oine  cm  intervals.  The  lowest  vertical  measurement  for  every  vertical 
profile  was  made  at  0.5  cm  from  the  bottom.  Since  the  probe  must  be  fully  submerged 
and  the  sampling  volume  is  53  mm  from  its  tip,  the  highest  measurement  attempted  was 
10  cm  from  the  free  surface. 

The  probe  was  put  in  a  tilted  position  for  all  measurements  in  front  of  the 
cylinder.  For  measurement  in  planes  perpendicular  to  the  flow  (i.e.  the  90°  plane  and  in 
the  wake  region)  the  probe  was  returned  to  the  vertical  position.  The  angle  of  the  tilt  was 
determined  before  conducting  the  experiment  so  that  the  probe  could  be  properly  located. 

Once  in  the  proper  position,  a  two  minute  record  was  acquired.  The  length  of  the 
record  was  a  compromise  between  taking  a  long  enough  sample  to  capture  the  longer 
period  flow  features,  such  as  the  wake  vortex  shedding,  getting  a  reliable  mean  flow 
vector  and  reducing  the  length  of  the  experiment  for  reasons  of  cost.  In  all, 
approximately  3900  records  were  acquired  totaling  over  130  hours  of  data. 

An  example  of  one  of  the  acquired  records  is  found  in  the  top  graph  of  Figure 
3.13.  Once  recorded  in  a  binary  file,  the  data  were  transferred  to  a  Pentium  computer 
where  it  was  reduced.  All  data  reduction  was  performed  using  the  Matlab  M-file 
dataredu.m.  This  program  is  listed  in  the  Appendix.  It  calls  the  subroutines  provided  by 
Sontek  that  read  the  binary  files  and  produce  ASCII  files  containing  the  desired 
information.  The  M-file  performs  the  following  calculations:  (1)  reads  the  velocity, 
position,  and  distance  to  the  boundary  from  the  binary  file,  (2)  rotates  the  velocity  data  to 
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flume  oriented  eoordmates  if  necessary,  (3)  removes  the  single  point  spikes  that  tend  to 
occur  in  the  raw  data,  (4)  removes  the  low  frequency  noise  associated  with  waves  in  the 
flume,  and  (5)  calculates  the  flow  statistics  including  the  three  mean  velocities,  the  terms 
of  the  Reynolds  stress  tensor  and  the  determinant  of  this  tensor. 

The  rotation  of  the  velocity  time  series  is  performed  by  examining  the  velocity 
vector  at  every  instant  of  time.  If  (ut,  Vt,  Wt)  is  the  velocity  vector  in  the  tilted  coordinate 
system  then  the  velocity  vector  in  the  flume  oriented  coordinate  system  is 
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0 
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where  Gxz  is  the  rotation  of  the  probe  in  the  X-Z  plane  and  6yz  is  the  rotation  of  the  probe 
in  the  Y-Z  plane.  A  comparison  of  the  raw  data  both  before  and  after  rotation  is 
presented  in  Figure  3.13. 

The  "smoothing"  of  the  data  involves  two  steps,  the  removing  of  the  single  point 
spikes  found  when  using  this  instrument  (see  the  middle  graph  in  Figure  3.13)  and  the 
removal  of  the  low  frequency  noise  of  the  free  surface  oscillation  that  occurs  in  both  the 
X  (along  flume)  and  Y  (across  flume)  directions.  The  time  record  is  scanned  for  data 
points  that  lie  three  standard  deviations  from  the  mean  value.  If  one  is  found,  the  data 
point  is  assigned  a  value  equal  to  the  mean  of  the  two  bracketing  values.  To  remove  the 
low  frequency  noise,  a  band  reject  filter  is  applied  to  the  Fourier  transform  of  the  velocity 
signal.  All  values  between  0  and  0.2  Hz  are  set  equal  to  zero.  The  inverse  transform 
recovers  the  original  signal  less  the  low  frequency  component.  Figure  3.13  shows  the 
final  result  after  filtering  and  rotating  the  original  raw  data. 
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Figure  3.13  Twenty  second  window  of  data  acquired  by  Sontek:  (top)  raw  data,  (center) 
raw  data  rotated  to  flume  oriented  coordinates,  and  (bottom)  rotated  and  filtered  data. 
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More  specific  information  concerning  the  manipulation  of  the  raw  data  is  found  in 
the  M-file  dataredu.m  in  the  Appendix.  The  time-averaged  velocity  measurements  from 
all  the  experiments  can  be  found  in  the  reference  Gosselin  and  Sheppard  (1997). 
3.2.3  Dye  Injection  Studies 

The  main  purpose  of  the  dye  injection  studies  was  to  gain  a  better  insight  into  the 
hydrodynamics  during  the  intermediate  stages  of  scour.  The  most  salient  features  of  the 
flow  field  are  the  horseshoe  vortex  and  the  wake  vortices.  As  such,  the  injections  were 
positioned  to  illuminate  these  mechanisms.  Photographs  of  the  horseshoe  vortex  for  the 
four  different  stages  of  scour  are  shown  in  Figure  3.14.  Photographs  of  the  wake  vortices 
are  shown  in  Figure  3.15. 

The  dye  studies  yielded  a  number  of  interesting  observations.  First,  the  horseshoe 
vortex  was  indeed  present  for  all  the  tests.  Figure  3. 14  clearly  shows  how  the  vortex  is 
initially  very  small  in  diameter.  As  the  hole  grows,  so  does  the  size  of  the  vortex. 
Accompanied  with  this  growth  is  an  apparent  weakening  of  the  strength  of  the  vortex. 
When  dye  was  injected  into  the  horseshoe  vortex  for  the  tests  with  the  1)  6  hours  and  2) 
24  hours  scour  holes,  the  vortex  appeared  to  climb  up  the  bottom  of  the  hole.  It  took  a 
number  of  seconds  before  the  dye  was  washed  away  from  this  region  as  opposed  to  tiie 
almost  instantaneous  removal  of  dye  for  the  flat  bed  and  T  =  1  hour  scour  hole.  This 
suggests  tiiat  the  velocities  next  to  the  bed  up  the  sides  of  the  hole  are  relatively  low  when 
compared  with  those  associated  with  the  scour  holes  at  earlier  stages  of  scour.  Figure 
3.14  (c)  and  (d)  show  Uie  dye  cUmbing  up  the  sides  of  the  holes  during  tiiis  investigation. 

The  dye  injection  studies  involving  the  wake  vortices  clearly  showed  the 
characteristic  von  Karman  vortex  sti-eet.  Unfortunately,  tiie  reflection  off  die  free  surface 
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made  photographing  tMs  phenomena  difficuU.  Two  of  the  best  are  shown  in  Figure  3.15. 
In  these,  flow  is  from  right  to  left.  The  curvature  of  the  streaks  and  the  concentration  of 
the  dye  gives  some  indication  of  these  vortices. 

3.3  Results  and  Discussion 

The  measurement  of  the  velocity  field  was  performed  in  vertical  planes.  The  time 
averaged  mean  velocity  is  presented  in  this  format.  The  ensuing  figures  show  the 
velocities  in  the  measurement  plane  in  the  form  of  velocity  vectors.  Velocities 
perpendicular  to  the  plane  are  represented  by  velocity  contours.  The  vectors  show  the 
direction  of  the  flow  in  that  plane  and  the  magnitude  of  the  velocity  is  proportional  to 
their  length.  The  length  of  the  vector  is  scaled  by  the  magnitude  of  the  largest  velocity 
vector  shown  in  the  figure.  This  magnitude  is  given  in  each  figure.  For  the  contours,  a 
positive  magnitude  indicates  flow  into  the  page. 

Measurements  are  presented  for  the  four  scour  cases.  For  each  figure,  the 
following  conventions  are  used:  (1)  the  coordinate  system  is  the  same  as  described 
previously,  (2)  the  origin  of  the  coordinate  system  is  at  the  center  of  the  cylinder  at  the 
initial  unscoured  bed  level,  (3)  planes  denoted  by  degrees  are  planes  rotated  from  the 
stagnation  plane  (Y  =  0  cm)  by  the  indicated  amount  passing  through  the  center  of  the 
cylinder,  and  (4)  planes  labeled  X  =  -b  and  X  =  -2b  are  Y-Z  planes  in  the  lee  of  the 
cylinder  at  distances  behind  the  cylinder  X  =  1  diameter  (=  -1 6.8  cm)  and  X  =  2  diameters 
(=  -33.7  cm).  For  the  initial  flat  bed  case,  measurements  were  made  in  the  planes  of  0°, 
45°  and  90°  from  tiie  stagnation  plane,  Y  =  1.8  cm,  Y  =  4.3  cm,  Y  =  7.3  cm,  X  =  -b  and  X 
=  -2b.  For  the  remaining  cases  measurements  were  made  in  the  planes  of  0°,  30°,  60°, 
and  90°  from  the  stagnation  plane  and  X  =  -b  and  X  =  -2b  (see  Figure  3.12). 
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Figure  3  . 14  Horseshoe  vortices  for:  (a)  t  =  0  hours, 
(b)  t  =  1  hour,  (c)  t  =  6  hours,  and  (d)  t  =  24  hours. 
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Figure  3.14  continued. 


88 


Figures.  15  (a)  and  (b)  Wake  vortices  (flow  is  from  right  to  left) 
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Figure  3. 16  Time  averaged  velocities  for  the  initial  flat  bed  case  in  the  planes:  (a)  0°,  (b) 
45°,  (c)  90°,  (d)  Y=1.8  cm,  (e)  Y=4.3  cm,  (0  Y=7.3  cm,  (g)  X=  -b,  and  (h)  X=-2b. 
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Figure  3.16  continued. 
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Figure  3.16  continued. 
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T  =  0  Hours,  X  =  -b 
Experimental  Data 


T  =  0  Hours,  X  =  -2b 
Experimental  Data 
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Figure  3.16  continued. 
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Figure  3.17  Time  averaged  velocities  for  the  T  =  1  hour  scour  hole  in  the 
planes:  (a)  0°,  (b)  30°,  (c)  60°,  (d)  90°,  (e)  X=-b,  and  (f)  X=-2b. 
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Figure  3.17  continued. 
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Figure  3.17  continued. 
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Figure  3.18  Time  averaged  velocities  for  the  T  =  6  hours  scour  hole  in 
planes:  (a)  0°,  (b)  30°,  (c)  60°,  (d)  90°,  (e)  X=-b,  and  (f)  X=-2b. 
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Figure  3.18  continued. 
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Figure  3.18  continued. 
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Figure  3.19  Time  averaged  velocities  for  the  T  =  24  hours  scour  hole  in 
the  planes:  (a)  0°,  (b)  30°,  (c)  60°,  (d)  90°,  (e)  X=-b,  and  (f)  X=-2b. 
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Figure  3.19  continued. 
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Figure  3.19  continued. 
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Included  in  each  of  the  figures  are  point  gauge  data  indicating  the  elevation  of 
bottom.  A  number  of  observations  can  be  made  from  the  preceding  figures.  Examining 
the  flat  bed  condition  in  Figure  3.16,  the  most  notable  feature  is  the  lack  of  a  strong 
horseshoe  vortex  at  the  leading  edge.  There  is  some  vorticity  at  the  leading  edge,  and 
perhaps  the  spatial  resolution  of  the  measurement  is  not  dense  enough  to  pick  up  a  more 
defined  vortex.  There  is  an  obvious  downflow  at  the  leading  edge  and  the  contours 
indicate  the  strong  redirection  of  the  flow  toward  the  side  of  the  cylinder.  The  90°  figure 
shows  the  substantial  acceleration  of  fluid  around  the  sides  of  the  cylinder  which  is  the 
most  influential  scouring  mechanism  during  this  phase. 

Figure  3.17  shows  the  velocity  vectors  after  one  hour  of  scouring.  During  this 
phase  of  scour,  the  hole  is  still  growing  quite  rapidly  (see  Figure  3.8).  More  fluid  has 
taken  a  downward  vertical  component  at  the  leading  edge  of  the  pile.  The  beginnings  of 
the  horseshoe  vortex  can  be  seen  at  the  very  bottom  of  the  hole.  There  is  no  separation  of 
the  flow  at  the  edge  of  the  hole.  This  may  be  the  reason  that  the  vortex  is  kept  small 
during  this  phase  of  scour. 

Figure  3.18  presents  the  flow  situation  for  the  scour  hole  after  6  hours  of  scour. 
The  scouring  process  at  this  point  has  slowed  significantly  and  has  started  a  gradual 
asymptotic  approach  toward  its  equilibrium  value.  Here,  there  is  a  clear  separation  at  the 
edge  of  the  scour  hole  and  a  well  defined  horseshoe  vortex,  The  vortex  remains  coherent 
as  it  is  swept  around  the  sides  of  the  cylinder.  The  flow  around  the  sides  from  the 
stagnation  plane  seems  strongest  in  the  center  of  the  vortex  inside  the  hole. 

Figure  3.19  shows  the  flow  field  as  the  hole  approaches  its  equilibrium  value  after 
24  hours  of  scour.  The  scouring  process  has  all  but  stopped  at  this  point.  The  velocity 
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vectors  show  a  very  coherent  vortex  in  the  stagnation  plane.  There  is  a  clear  separation  at 
the  upstream  edge  of  the  scour  hole.  The  vortex  can  be  seen  narrowing  as  it  is  pulled 
around  the  edge  of  the  cylinder.  It  also  is  drawn  away  from  the  cylinder  and  up  the  side 
of  the  scour  hole  as  it  moves  around.  The  intensification  of  the  vortex  is  difficult  to  judge 
since  the  planes  are  shown  at  different  angles  to  the  main  flow  direction. 

The  sequence  of  figures  can  be  summarized  as  follows:  initially,  scour  is 
performed  primarily  by  the  acceleration  of  the  fluid  around  the  sides  of  the  cylinder  and 
the  downflow  at  the  leading  edge.  As  the  scour  hole  deepens,  more  fluid  attains  a 
downward  component  at  the  leading  edge.  This  fluid  impinges  on  the  bed  and  is  swept 
around  the  side  of  the  cyhnder.  After  1  hour,  a  well  defined  horseshoe  vortex  has  not  yet 
appeared.  Upon  further  deepening  of  the  hole,  the  fluid  eventually  separates  at  the 
upstream  edge  of  the  hole.  This  allows  the  creation  of  a  coherent  horseshoe  vortex.  The 
vortex  is  swept  around  the  sides  of  the  cylinder  were  it  narrows  in  cross-section  and 
climbs  up  the  edges  of  the  hole.  Continued  deepening  of  the  hole  is  accompanied  by  the 
expansion  of  the  horseshoe  vortex  to  fill  the  hole  and  a  more  coherent  vortex  structure 
around  the  edges  of  the  pile.  The  horseshoe  vortex  is  easily  recognized  at  this  point  due 
in  part  to  the  low  velocities  found  inside  the  hole  as  the  equilibrium  condition  is  neared. 

Examination  of  the  wake  region  shows  interesting  behavior  for  the  time  averaged 
data.  Looking  at  the  figures  for  the  plane  X  =  -b,  there  is  a  well  defined  wake  region  of 
low  velocity  magnitude  in  the  lee  of  the  cylinder  for  the  flat  bed  condition.  There  is 
relatively  small  flow  laterally  for  this  case  with  a  maximum  magnitude  of  7.4  cm/p. 
Once  scour  starts,  however,  this  maximum  doubles  in  size.  This  may  be  due  to  the  fact 
that  the  scour  hole  itself  redirects  flow  toward  the  lee  of  the  cylinder.  As  the  hole 
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expands,  the  lower  velocities  found  in  the  expanded  cross-section  translate  to  a  lower 
magnitude  of  the  velocity  component  in  the  measurement  plane.  For  the  planes  X  =  -2b, 
the  flow  situation  is  more  similar  in  structure  between  different  times  with  similar 
differences  in  flow  magnitudes.  The  upwelling  associated  with  the  wake  vortices  can  be 
seen  in  all  cases. 

It  is  also  instructive  to  view  the  data  in  horizontal  planes:  Figures  3.20  through 
3.22.  Here  the  figures  are  organized  by  plane  rather  than  by  time.  Three  situations  are 
shown.  Figure  3.20  presents  the  bottom  velocities.  Although  an  attempt  was  made  to 
obtain  measurements  5  mm  from  the  bottom  for  every  vertical  profile,  sometimes 
logistics  and  the  fact  that  the  ADV  does  not  work  well  near  boundaries  prohibited  this. 
As  a  result,  the  figures  show  the  data  closest  to  the  bottom  for  each  vertical  profile.  The 
fact  that  the  distance  to  the  bottom  varies  for  each  measurement  should  be  kept  in  mind 
when  viewing  these  figures.  Figures  3.21  and  3.22  show  the  velocities  in  the  plane  Z  =  1 
and  5  cm. 

The  bottom  velocity  vectors  illustrate  the  ti^ansition  from  the  scour  due  to 
accelerated  flow  around  tiie  sides  of  the  cylinder  to  the  scour  due  to  the  horseshoe  vortex. 
Figure  3.19  (a)  and  (b)  show  the  strong  acceleration  of  tiie  fluid  near  die  pile  reaching  a 
maximum  velocity  at  tiie  90°  plane.  In  Figures  3.19  (c)  and  (d)  tiie  flow  in  tiie  0°  and  30° 
plane  show  the  flow  up  tiie  side  of  tiie  hole  in  the  opposite  direction  to  the  upsti-eam 
velocity.  This  flow  is  the  bottom  velocity  associated  with  the  horseshoe  vortex.  Figures 
3.21  and  3.22  also  show  the  two  dimensional  acceleration  of  the  fluid  around  the  edges  of 
the  cyhnder  as  well  as  the  retardation  of  the  flow  in  the  wake  region. 
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Figure  3.20  Velocities  at  the  bed  at  times  :  (a)  T  =  0  hours, 
(b)  T  =  1  hour,  (c)  T  =  6  hour,  and  (d)  T  =  24  hours. 
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Figure  3.20  continued. 
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Figure  3.21  Velocities  for  the  plane  Z  =  +1  cm  at  times:  (a)  T  =  0 
hours,  (b)  T  =  1  hour,  (c)  T  =  6  hour,  and  (d)  T  =  24  hours. 
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T  =  6  Hours,  Z  =  1  cm 


Largest  vector  =  31 .9cm/s 
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Figure  3.21  continued. 
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T  =  0  Hours,  Z  =  5  cm 
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Figure  3.22  Velocities  for  the  plane  Z  =  +5  cm  at  times:  (a)  T  =  0 
hours,  (b)  T  =  1  hour,  (c)  T  =  6  hour,  and  (d)  T  =  24  hours. 
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T  =  6  Hours,  Z  =  5  cm 


Largest  vector  =  32.83cm/s 
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Figure  3.22  continued. 


Ill 


An  attempt  was  also  made  to  examine  the  frequency  of  the  wake  vortex  shedding. 
The  region  examined  was  directly  behind  the  left  edge  of  the  cylinder,  one  diameter  back 
from  the  center  of  the  cylinder  and  high  in  the  water  column.  The  Fourier  transform  of 
the  velocity  time  signal  at  this  point  was  analyzed  for  each  of  the  four  cases.  The  peak  in 
the  plot  of  the  magnitude  of  the  Fourier  transform  should  correspond  to  the  frequency  of 
the  vortex  shedding.  According  to  Schlichting  (1955),  the  Strouhal  number 
corresponding  to  a  Reynolds  number  of  4.1  x  10"*  (the  Reynolds  number  associated  with 
the  cy Under  for  these  experiments)  is  S  =  0.19.  This  translates  to  a  shedding  frequency  of 
0.28  Hz  for  the  conditions  of  this  test.  For  the  flat  bed  case,  the  plot  of  the  magnitude  of 
the  Fourier  transform  of  the  velocity  in  the  X-direction  is  given  in  Figure  3.23. 
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Figure  3.23  Magnitude  of  the  Fourier  transform  of  the 
velocity  in  the  X-direction  in  the  wake  region. 
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The  value  of  the  frequency  at  the  peak  is  0.2808  Hz.  For  the  cases  of  T  =  1 ,  6, 
and  24  hours,  the  values  of  the  frequency  at  the  peak  are  0.2441  Hz,  0.2319  Hz,  and 
0.2197  Hz,  respectively.  At  24  hours,  there  is  a  22%  reduction  in  the  shedding  frequency. 
As  Melville  (1975)  pointed  out,  this  reduction  is  due  to  the  decreased  velocities  that  the 
pile  experiences.  The  increase  in  the  scour  depth  leads  to  a  larger  cross-sectional  area  and 
hence  a  lower  velocity.  Since  the  shedding  frequency  is  directly  proportional  to  the 
velocity  and  the  Strouhal  number  is  relatively  constant  for  this  range  of  Reynolds 
numbers,  the  frequency  of  shedding  decreases. 


CHAPTER  4 
NUMERICAL  EXPERIMENTS 


This  chapter  describes  the  computer  simulations  performed  as  part  of  this 
investigation.  It  begins  with  a  description  of  the  software,  FL0W-3D®,  used  for  the 
simulations.  Next,  a  description  of  the  specific  features  used  for  the  computational 
experiments  performed  is  presented.  This  chapter  will  conclude  with  a  presentation  and 
discussion  of  the  specific  experiments  and  their  results. 

4.1  Computer  Software 

FL0W-3D®  is  a  commercially  available  software  package  developed  by  Flow 
Science,  Inc.  The  software  is  capable  of  analyzing  a  range  of  fluid  dynamic  and  thermal 
phenomena.  It  is  equipped  to  model  three-dimensional,  time  dependent  flows  as  well  as 
handle  both  complex  geometric  forms  and  free  surfaces.  It  also  provides  a  number  of 
options  for  turbulent  closure.  The  software  is  also  available  for  a  number  of  platforms 
ranging  from  Pentium  PCs  to  super  computers. 
4.1.1  Program  Components 

FL0W-3D®  consists  of  four  individual  programs.  The  preprocessor 
(PREP3D.EXE)  simplifies  the  setup  of  the  problem  by  translating  simple  input 
statements  into  the  comprehensive  numerical  description  of  the  desired  flow  situation  to 
be  read  by  the  main  processor.  The  main  processor  (HYD3D.EXE)  reads  the  input  files 
generated  by  the  preprocessor  and  produces  the  numerical  solution  to  the  problem.  This 
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program  incorporates  not  only  the  algorithms  for  producing  solutions  to  many  flow 
situations,  but  also  contains  the  abiUty  to  optimize  the  solution  times  through  the  choice 
of  convergence  criteria  and  time  step  sizes.  The  post-processor  (FLSCON.EXE)  reads 
the  output  files  created  by  the  main  processor  and  manipulates  the  data  in  the  manner 
desired  by  the  user  (i.e.  images,  ASCII  file  output,  etc.).  Finally,  a  plotting  package 
(PLTFSI.EXE)  is  included  that  reads  the  image  files  created  by  both  the  pre  and  post- 
processors and  displays  them. 

Execution  of  a  numerical  experiment  using  FL0W-3D®  begins  with  the 
preparation  of  m  input  file  called  prepin.inp.  The  input  file  consists  of  a  title  block 
followed  by  a  list  of  input  blocks  associated  with  some  aspect  of  the  code  operation  or  the 
description  of  the  problem.  These  blocks  are  lists  of  variables  organized  by  which  aspect 
of  the  code  they  address:  for  example,  initial  conditions,  boundary  conditions,  fluid 
properties,  etc.  An  example  of  the  format  of  a  typical  input  file  used  in  this  investigation 
is  found  in  Figure  4.1.  This  format  allows  for  quick  setup  of  the  simulation  as  well  as 
fast  modification  of  existing  input  files. 

The  input  file  is  read  by  the  preprocessor.  The  preprocessor  converts  the  data  into 
a  form  suitable  for  the  main  processor.  As  dictated  by  the  input  file,  it  initializes  the  fluid 
properties,  sets  the  initial  conditions,  fixes  the  boundary  conditions,  generates  the 
obstacles  and  defines  all  the  computational  control  parameters.  It  uses  an  automatic  mesh 
generator  for  rapid  creation  of  a  computational  grid  given  a  number  of  specifications  by 
the  user.  Rather  than  using  a  body  fitted  coordinate  system,  the  software  uses  an 
orthogonal  grid.  Objects  are  handled  using  the  FAVOR  (fractional  area/volume)  method 
that  will  be  described  in  the  next  section.  Upon  completion,  the  preprocessor  creates  a 
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number  of  input  files  for  the  main  processor.  Included  among  them  is  a  file  containing  a 
series  of  images  that  can  be  read  by  the  plotting  program.  This  permits  quick  and  easy 
debugging  of  flow  experiments. 


Flow  around  a  Cylinder  at  Equilibrium  Scour  Depth 
cylinder  diameter  =  0.1 68m  =  6  5/8  in 
water  depth  =  34.83  in 


&xput 

remark= 'units  are  cgs', 
dtmin=1.0e-12,  twfm=20., 

/ 


ipdis=l. 


&limits 
irpr=2, 

/ 

&props 

rhof=l., 

/ 

&bcdata 
wl=6, 

/ 

&mesh 

nxcelt=40, 

/ 

&obs 

nobs=2, 

/ 

&fl 

flht=34.83, 

/ 

&temp 

tempi=293., 

/ 

&grafie 
nvplts=3, 

/ 


jbkpr=2,  ktpr=2, 
mui=0.01, 

wr=3, 


wf=l 


px(l)=-40.0,  px(2)=-3 1.748, 
iofo(l,l)=l,  zl(l)=-14.. 


contpv(l)='p',  xvl(l)=-25., 


Figure  4. 1  Example  of  the  format  of  a  prepin.inp  file. 
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The  next  step  in  the  execution  of  the  software  is  performed  by  the  main  processor. 
This  program  calculates  the  actual  solution  to  the  problem.  It  uses  a  three-dimensional, 
finite  difference  method  to  calculate  the  time-dependent  solution.  The  equations  solved 
are  the  conservation  laws  {qt  mass,  momentum  and  energy.  The  difference  equations  use 
a  fixed  Eulerian  grid  of  non-uniform  rectangular  control  volumes.  A  fractional  volume  of 
fluid  function  (VOF)  is  used  to  handle  free  surfaces  and  interfaces.  This  function  will  be 
discussed  in  the  next  section.  Upon  completing  the  solution,  the  main  processor  creates  a 
number  of  output  files.  From  these  output  files  the  user  has  the  ability  to  restart  a 
simulation  from  the  end  of  a  previous  simulation. 

The  post-processor  uses  the  data  files  from  the  main  processor.  It  has  the  ability 
to  create  a  number  of  graphical  presentations  of  the  data  as  well  as  output  the  data  in 
ASCII  form.  Graphical  capabihties  include  vector  plots,  contour  plots,  time  history  plots, 
contour  surface  plots,  fluid  surface  plots  and  marker  particle  distributions.  The  post- 
processor creates  a  file  containing  the  requested  graphics  to  be  read  by  ttie  plotting 
program. 

The  image  display  program  reads  the  binary  files  created  by  both  the  pre  and  post- 
processors. It  creates  hardware- specific  displays  of  the  data  and  is  primarily  used  for 
graphical  display  data  on  the  computer  monitor.  It  has  the  option  of  creating  output  files 
of  the  individual  displays  in  both  PostScript  and  bitmap  formats. 
4.1.2  Program  Features 

The  differential  equations  solved  by  the  main  processor  are  written  in  terms  of 
Cartesian  coordinates.  The  solution  mesh  consists  of  arrays  of  non-uniform  box  shaped 
elements.  FL0W-3D®  also  has  the  capacity  of  computation  in  a  cylindrical  coordinate 
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system.  For  this  case,  terms  are  added  to  the  equations  of  motion.  In  this  discussion, 
only  the  Cartesian  formulations  will  be  considered.  For  the  finite  difference  analysis,  the 
following  conventions  apply:  scalar  fluid  properties,  such  as  fluid  fraction,  pressure, 
density,  viscosity,  etc.,  are  evaluated  at  the  center  of  the  control  volume  and  velocities  are 
evaluated  at  the  faces  between  adjoining  control  volumes  in  the  X  direction  for  the  u 
velocity,  in  the  Y  direction  for  the  v  velocity  and  in  the  Z  direction  for  the  w  velocity. 
The  convention  is  illustrated  in  Figure  4.2. 


In  Figure  4.2,  F  is  the  fluid  fraction,  q  is  the  turbulent  kinetic  energy,  D  is  the  turbulent 
dissipation,  Vf  is  the  fraction  of  the  volume  available  for  fluid  flow  and  Ax,  Ay,  and  Az  are 
the  fractional  areas  of  the  X,  Y,  and  Z  faces. 

Before  presentation  of  the  equations  of  motion,  it  is  helpful  to  review  how  the 
software  handles  obstacles  and  complex  geomefric  regions.  The  code  makes  use  of  a 
formulation  known  as  the  FAVOR  method  (Fractional  AreaA'^olume  Obstacle 
Representation).  This  technique  gives  a  simple,  accurate  way  to  avoid  "saw-tooth" 


Figure  4.2  Location  of  variables  in  a  mesh  cell. 
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boundaries  and  gives  a  good  geometric  description  of  obstacles  even  in  relatively  coarse 
grids.  Once  the  user  has  input  the  desired  obstacle,  the  code  determines  where  the  object 
is  embedded  in  the  computational  mesh.  For  each  control  volume,  it  calculates  what 
fraction  of  the  volume  the  object  occludes  as  well  as  what  fractional  area  on  each  face  of 
the  volume  is  blocked.  In  this  manner,  the  code  decides  what  faces  of  the  cell  are  open  to 
the  flow  and  how  much  fluid  is  contained  in  that  cell.  To  illustrate  this  point,  Figure  4.2 
contains  a  cell  with  an  obstacle  intersecting  it.  Viewing  the  right  face  of  the  cell  from  the 
positive  Y-direction  yields  Figure  4.3. 

Ay  (i,j,k) 


A,(i,j,k-1) 

Figure  4.3  Illustration  of  FAVOR  method. 
Knowing  that  the  fraction  of  the  area  of  the  illustrated  (Y)  right  face  is  blocked  and  the 
fact  that  the  front  (X)  face  of  the  cell  directly  behind  is  blocked  situates  the  obstruction  at 
the  back  right  of  the  cell.  Add  to  this  the  fact  that  the  top  (Z)  face  of  the  cell  directly 
below  is  blocked  situates  the  obstruction  at  the  lower,  back,  right  comer.  These  fractional 
areas  and  volumes  are  included  in  the  finite  difference  equations. 
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The  equations  solved  by  the  main  processor  are  the  conservation  laws  for  mass, 
momentum  and  energy.  The  mass  conservation  law  has  the  following  form  for 
incompressible  flow: 


a(uAj    3(vaJ    a(wAj  RSOR 


dx  dy  dz  p  ' 

where  RSOR  is  a  mass  source  term  and  u,  v,  and  w  are  the  x,  y,  and  z  components  of 

velocity  respectively.  The  momentum  conservation  equations  are  summarized  in 

Equation  4.2. 


au   1  f     au      au       au  1     i  ap  rsor 

— +— "^uA  — +  vA  — +  wA  — ^  =  — — +  G,  +f^-b^  — u 

at  Vf  [  ^  ax      ay       azj    p  ax    ^   "    "  pv 

av     1  f       dv  dy  dv]       1  ap  RSOR 

iuA^  — +  vA^  — +  wA^  — ^  =  — — +  G  +f  -b  - — — V     .  (4.2) 
at    Vf  [     "  ax        ^  ay        'azJ       pay      y     y     y      pv  ^  ^ 

aw   1  f    aw      aw       awl    i  ap  rsor 

— +— <^uA,— +vA  — +  wA,  — ^  =  — — +  G,  +f,  -b,  - — — w 

at  Vf  [     ax     ^  ay       az  j    p  az    ^   ^    ^  py 


(Gx,  Gy,  Gz)  are  body  accelerations,  (fx,  fy,  fz)  are  viscous  accelerations,  (bx,  by,  bz)  are 
flow  losses  across  porous  baffles  and  the  final  terms  account  for  the  injection  of  mass. 

The  software  gives  a  number  of  options  for  turbulence  closure.  These  options 
include  the  Prandtl  mixing  length,  a  turbulent  energy  model,  a  two  equation  k-e  model, 
the  Renormalized  Group  Theory  (RNG)  model,  and  a  large  eddy  simulation  model.  The 
choice  for  turbulence  closure  for  this  investigation  will  be  discussed  later. 

The  code  handles  free  surfaces  using  a  method  known  as  the  Volume  of  Fluid 
(VOF)  technique  pioneered  by  Hirt  and  Nichols  (1981).  This  technique  consists  of  three 
components:  a  method  for  finding  the  free  surface,  an  algorithm  for  tracking  the  free 
surface  as  a  sharp  interface  moving  through  the  computational  mesh  and  a  process  for 
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applying  boundary  conditions  to  the  surface.  The  method  makes  use  of  the  simple 
principle  of  assigning  a  single  variable  F  (fluid  fraction)  to  each  cell  that  has  a  value  of 
1.0  if  the  cell  is  occupied  by  fluid  and  a  value  of  0.0  if  the  cell  is  completely  empty. 
Therefore,  if  the  cell  has  a  value  of  F  between  0  and  1  then  the  cell  contains  a  free 
surface.  In  addition,  the  normal  to  the  surface  can  be  calculated  from  the  direction  in 
which  F  changes  most  rapidly.  The  evolution  of  the  F  field  is  dictated  by  the  differential 
equation  below. 

aF    aF    ap  aF 

This  equation  dictates  that  F  moves  with  the  fluid.  In  this  manner,  FL0W-3D®  offers  a 
method  for  tracking  the  movement  of  the  free  surface  with  a  minimum  of  storage 
requirements. 

Boundary  conditions  on  solid  surfaces  are  handled  individually  on  each  surface  if 
multiple  surfaces  are  found  within  a  control  volume.  First  the  cell  is  examined  to 
determine  whether  a  face  has  a  fractional  area  less  than  unity.  If  so,  the  acceleration  due 
to  wall  shear  is  calculated.  For  example,  an  X  directed  velocity  bounded  on  the  bottom 

au 

would  have  a  shear  stress  approximated  by  the  code  using  the  term  A|li— .  The  variable 

dz 

A  is  included  to  account  for  turbulent  boundary  stresses.  For  laminar  flow  A  has  a  value 
of  1.0.  If  the  flow  is  turbulent,  a  "law  of  the  wall"  profile  exists  and  the  factor  A  must 
depend  upon  the  local  Reynolds  number.  FL0W-3D®  uses  a  1/7*  power  law 
approximation  to  the  logarithmic  profile.  The  factor  A  is  then  described  by  the  following 
equation: 
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1.0        if    Rl  <  157.44 
0.0225Rl'    if    Rl>  157.44  ^"^^ 

where  Rl  is  the  local  Reynolds  number  based  on  the  size  of  the  cell  in  the  direction  of  the 
gradient.  The  use  of  the  1/7"^  power  law  rather  than  a  logarithmic  profile  eliminates  the 
need  to  solve  a  transcendental  equation. 

Mesh  boundaries  can  take  one  of  six  different  conditions:  symmetry  boundaries, 
wall  boundaries,  continuativc  boundaries,  periodic  boundaries,  fixed  pressure  boundaries 
and  fixed  velocity  boundaries.  Symmetry  boundaries  are  used  across  planes  where  the 
flow  parameters  are  mirror  images.  As  such,  the  normal  velocities  are  zero  across  these 
planes.  Wall  boundaries  also  have  zero  normal  velocity  and  have  the  option  of  applying  a 
no  slip  condition  on  them.  Continuative  boundaries  are  not  physical  boundaries  per  se, 
but  rather  numerical  approximations  to  a  free  outflow  boundary.  Periodic  boundaries 
specify  that  the  flow  conditions  are  equal  at  two  opposing  mesh  faces.  They  therefore 
must  be  applied  in  pairs.  Fixed  pressure  and  fixed  velocity  boundaries  impose  a  specified 
pressure  or  velocity  over  the  entire  mesh  boundaries.  The  user  has  the  option  of  making 
these  conditions  time  dependent.  Unfortunately,  in  the  software's  commercial  form  there 
is  no  option  for  applying  a  velocity  boundary  condition  that  varies  spatially  over  the  mesh 
boundary  (i.e.  a  velocity  profile).  This  seriously  hampers  its  ability  to  simulate  open 
channel  flows. 

Finally,  obstacle  generation  using  FL0W-3D®  is  one  of  its  most  convenient 
features.  Since  the  obstacles  are  independent  of  the  computational  mesh,  their  definition 
may  be  made  in  terms  of  Cartesian  coordinates  rather  than  in  its  relationship  to  the  grid. 
As  a  result,  the  software  lets  the  user  define  obstacles  in  terms  of  functions  of  X,  Y,  and 
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Z.  The  software  also  contains  a  number  of  "canned"  objects  such  as  blocks,  spheres  and 
cylinders.  In  addition,  complex  geometric  shapes  may  be  configured  by  overlapping  a 
number  of  individual  objects.  Objects  are  defined  by  whether  they  are  solid  or  a  hole  in 
an  existing  object.  For  example,  a  cylindrical  hole  in  a  sphere  is  constructed  by  first 
defining  a  sphere  and  then  defining  a  cylinder  intersecting  the  sphere  and  telling  the 
software  that  the  cylinder  is  a  hole.  The  ease  of  defining  obstacles  makes  the  chore  of 
constructing  a  complex  shape,  such  as  a  scour  hole  around  a  cylindrical  pile,  fast  and 
efficient. 

4.2  Computational  Simulations 

The  purpose  of  performing  computational  experiments  is  twofold:  first,  to  test  the 
ability  of  the  software  to  accurately  model  the  flow  in  a  scour  hole  by  comparing  the 
results  of  the  simulations  with  the  experimental  measurements,  and  second,  if  the  first 
criteria  is  met,  to  use  the  program  to  evaluate  the  flow  patterns  in  the  scour  hole  as  it 
evolves  with  time.  This  section  is  divided  into  three  parts.  The  first  section  addresses  the 
individual  program  options  employed  for  the  simulations  including  boundary  conditions 
and  turbulence  closure.  The  second  describes  the  simulations  that  model  the  fixed  bed 
physical  experiments.  The  third  section  outlines  the  simulations  performed  to 
characterize  the  time  dependent  evolution  of  the  scour  hole.  All  the  computational 
experiments  in  this  investigation  were  performed  using  the  FL0W-3D®  software  running 
on  a  90  MHz  Pentium  PC  with  32  Mbytes  of  RAM. 
4.2.1  Program  Options 

This  section  will  explain  the  choices  of  the  many  options  offered  by  the  software, 
used  in  this  work.  Since  the  purpose  of  the  simulations  was  to  examine  the  mechanisms 
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of  scour,  computationally,  care  was  taken  to  focus  on  the  main  scouring  agents,  the 
downflow  at  the  front  of  the  cylinder  and  the  horseshoe  vortex.  As  such,  resolution  of 
some  of  the  secondary  flows,  such  as  the  wake  region  and  the  surface  roller,  was  of  lesser 
importance.  Some  of  the  choices  of  modeling  parameters  reflect  this  hierarchy.  The 
organization  of  this  section  will  follow  that  of  the  input  file  hsted  in  Figure  4.1  (i.e.  the 
section  headings  beginning  with  the  ampersand,  &).  First,  choices  for  the  physical  and 
computational  parameters  are  discussed  (&xput  section).  Next,  the  boundary  condition 
parameters  are  reviewed  (&bcdata  section).  Finally,  the  computational  mesh  and  the 
creation  of  the  obstacles  are  treated  (&mesh  and  &obs  sections).  For  all  the  simulations 
performed,  the  fluid  employed  had  a  density  of  1.0  g/cm^  and  a  dynamic  viscosity  of  0.01 
g/(cm  s). 

The  physical  situation  modeled  by  the  simulations  is  horizontal  flow  over  a  flat 
bed.  Gravity  acts  vertically  downward  at  980.6  cm/sec^  creating  a  vertical  hydrostatic 
pressure  distribution.  The  option  for  wall  shear  stress  is  used  so  that  there  is  a  no-slip 
condition  on  all  boundaries  designated  as  walls  and  all  obstacles.  The  choice  for  the  time 
Of  simulations  was  20  seconds.  Given  a  "hot-start"  to  the  model,  20  seconds  allowed 
sufficient  time  for  the  model  to  reach  an  equihbrium  condition  at  least  in  the  regions  of 
interest.  The  "hot-start"  will  be  discussed  later  in  the  boundary  conditions  section. 

Turbulence  closure  was  achieved  using  the  RNG  model.  The  Renormalized 
Group  Theory  model  for  turbulence  was  developed  by  Yakhot  and  Orzag  (1986)  and  is 
considered  to  be  an  improvement  on  the  well  established  two  equation  k-e  model.  In 
general,  mixing  length  and  eddy  viscosity  models  require  a  high  degree  of  empiricism  in 
order  to  model  phenomena  such  as  near  wall  turbulence  damping  and  boundary  layer 
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transition.  The  RNG  approach  applies  a  formal  analysis  to  turbulence  that  systematically 
removes  the  small  scales  of  turbulence  up  to  a  resolvable  scale  and  then  includes  the 
effect  of  these  scales  through  an  effective  or  renormalized  viscosity.  The  result  is  a 
turbulence  model  that  includes  a  minimum  of  experimentally  adjustable  parameters.  This 
method  shows  improved  results  in  modeling  swirling  flows,  flows  in  regions  of  high 
distortion,  wakes  and  separated  flows.  In  addition  to  these  benefits,  the  RNG  model  also 
provided  smaller  computation  times  than  the  less  complex  k-e  and  Prandtl  mixing  length 
models  in  some  instances.  For  the  simulation  of  the  scour  hole  after  24  hours,  the  RNG 
model  completed  20  simulated  seconds  in  13.48  hours  while  the  k-e  and  Prandtl  models 
took  16.01  and  14.72  hours  respectively.  In  comparison  of  the  results  from  the  different 
models,  it  was  found  that  the  results  of  the  mean  flow  quantities  were  largely  similar 
except  in  the  wake  region  where  the  RNG  model  provided  slightly  improved  results  when 
compared  with  the  measured  data.  Hence,  since  the  RNG  model  rehed  on  fewer 
empirical  parameters  and  actually  provided  an  improvement  over  the  k-e  model  in 
computational  time,  it  was  employed  for  the  turbulence  closure  scheme  for  the 
simulations  in  this  investigation. 

The  boundary  conditions  on  the  computational  domain  were  chosen  to  minimize 
the  number  of  grid  points  while  preserving  the  physics  of  the  model  To  begm,  the  left 
side  of  the  domain  was  designated  as  the  fixed  velocity  boundary.  Therefore,  flow 
through  the  mesh  occurs  from  left  to  right.  It  was  mentioned  earlier  that  the  software  did 
not  allow  for  a  spatially  varying  boundary  condition.  Given  this  constraint,  the  only  way 
to  produce  a  boundary-layer,  open-channel  flow  was  to  specify  a  constant  velocity  across 
the  boundary  equal  to  the  desired  mean  value  and  then  extend  the  grid  sufficiently  far 
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upstream  to  allow  sufficient  distance  for  the  boundary  layer  to  develop.  It  should  be 
noted  at  this  point  that  there  is  no  computational  method  for  giving  comphcated 
boundaries  or  obstacles  roughnesses  for  turbulence  generation  purposes.  The  obvious 
drawback  to  this  method  is  that  it  causes  a  dramatic  increase  in  the  number  of 
computational  cells  required  and  the  upstream  velocity  profile  that  is  produced  is  for  a 
hydraulicly  smooth,  flat  plate  rather  than  one  with  the  desired  roughness. 

To  resolve  this  inadequacy,  Flow  Science,  Inc.,  was  contacted  directly  and  they 
agreed  to  customize  the  code  for  this  investigation  to  allow  for  the  specification  of  a 
vertically  varying  velocity  profile  that  mimics  the  behavior  of  the  experimental  data.  It 
was  decided  that  the  best  way  to  do  this  was  in  a  piece-wise  fashion:  the  near-wall 
behavior  is  represented  by  a  power  function  (u  =  Az^)  and  the  upper  region  of  the  flow  is 
represented  by  a  straight  line  (u  =  Cz  -I-  D).  The  constants  A  through  D  and  the  point  of 
transition  from  one  function  to  the  other  are  specified  by  the  user.  Figure  4.4  shows  the 
velocity  profile  used  for  the  simulation  of  the  flat  bed  (l  ^  0  hours)  condition. 

The  circles  show  the  experimental  data  at  the  same  X-position  were  the  boundary 
condition  is  applied.  For  the  simulations  that  model  the  physical  experiments,  the 
constants  for  the  functions  and  the  point  of  transition  were  found  using  a  multivariable 
least  squares  fit  to  the  appropriate  experimental  data.  The  implementation  of  the  user 
specified  boundary  condition  significantly  reduces  the  number  of  grid  points  used  since 
the  boundary-layer  velocity  profile  is  "hardwired"  into  the  boundary  condition.  Although 
there  is  still  no  way  to  introduce  bed  roughness  into  the  code,  the  inclusion  of  the  velocity 
profile  reduces  the  influence  of  this  discrepancy  on  the  mean  velocities  away  from  the 
boundaries. 
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Figure  4.4  Best  fit  user  specified  boundary  condition  to  experimental 
data  for  the  (lat  bed  (t  -  0  hours)  computer  simulation. 


The  downstream  (right)  boundary  condition  is  a  continuative  boundary  that  allows 
the  flow  to  leave  the  computational  mesh  (see  Figure  4.5).  The  bottom  boundary  is  a  wall 
boundary.  The  front  boundary  is  a  symmetry  condition.  This  stipulates  that  there  is  no 
flow  through  the  front  boundary.  It  also  models  the  flow  situation  where  there  is  an  array 
of  cylinders  perpendicular  to  the  flow  at  a  spacing  equal  to  the  width  of  the  grid.  Care 
was  taken  to  make  the  mesh  wide  enough  to  eliminate  the  influence  of  this  boundary 
condition.  The  back  boundary  was  also  made  a  symmetry  boundary  condition  where  the 
plane  of  symmetry  cuts  through  the  center  of  the  cylinder  and  scour  hole.  In  other  words, 
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the  plane  of  symmetry  is  coincident  with  the  stagnation  plane  of  the  cylinder.  The  main 
reason  for  doing  this  was  to  cut  by  Vi  the  required  number  of  cells.  The  obvious 
drawback  to  this  is  that  it  is  no  longer  possible  to  model  the  wake  region  adequately.  It  is 
not  possible  to  capture  alternately  shedding  vortices  by  using  a  symmetry  condition 
through  the  center  of  the  cyhnder.  Imposing  this  condition  did  not  affect  the 
hydrodynamics  at  the  leading  edge,  however.  As  a  result,  it  was  decided  that  the  savings 
accrued  in  computation  time  justified  the  imposition  of  this  condition. 

The  final  boundary  condition  is  the  fixed  pressure  condition  at  the  top  of  the  mesh 
boundary.  The  pressure  at  this  boundary  is  atmospheric.  An  additional  condition  was  put 
on  this  face.  Since  this  face  is  entirely  out  of  the  water,  it  is  possible  to  impose  the 
condition  that  the  fluid  fraction  at  the  top  face  be  equal  to  zero.  The  reason  for  doing  this 
is  to  speed  up  computation  time  for  the  "hot-start."  The  "hot-start"  refers  to  starting  the 
simulation  by  giving  all  the  fluid  in  the  computational  domain  a  velocity  equal  to  the 
average  velocity  of  the  boundary  condition  on  the  left  side  of  the  mesh.  This  eliminates 
having  to  ramp  the  velocity  up  from  zero  to  the  desired  value.  It  also  significantly 
reduces  the  computational  time.  One  problem  associated  with  this  technique,  however,  is 
that  the  fluid  located  next  to  the  cylinder  is  given  a  relatively  large  velocity  and  has 
nowhere  to  go  but  straight  up.  The  splash  that  occurs  at  the  start  of  the  simulation  is  very 
computationally  intensive.  By  imposing  the  fluid  fraction  condition  on  the  top  of  the 
mesh,  any  fluid  that  splashes  high  enough  to  reach  the  top  is  simply  removed  from  the 
computational  domain.  This  helps  speed  up  the  simulation.  The  use  of  20  seconds  for  a 
simulation  time  was  chosen  by  trial  and  error  as  the  point  at  which  the  flow  achieved 
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equilibrium  after  the  "hot-start."  A  review  of  the  boundary  conditions  for  all  the 
simulations  run  is  contained  in  Figure  4.5. 

The  meshes  utihzed  in  the  simulations  varied  sUghtly  from  experiment  to 
experiment  due  to  the  fact  that  mesh  generation  is  controlled  automatically  by  the 
software.  In  creating  the  mesh,  the  user  specifies  the  total  number  of  nodes  and  the 
distance  between  the  opposing  boundaries  for  each  of  the  three  directions.  In  addition, 
the  user  is  allowed  the  option  of  specifying  the  placement  of  a  solution  plane  between  the 
boundaries,  the  spacing  on  either  side  of  this  plane,  and  the  number  of  solution  planes 
contained  between  specified  planes.  To  provide  some  consistency  between  numerical 
experiments,  a  few  rules  were  followed  in  the  mesh  section  of  the  prepin.inp  file.  First, 
the  center  of  the  cylinder  was  located  at  (X,Y)  =  (0,0),  and  the  flow  is  from  left  to  right  in 
the  positive  X  direction.  Since  the  radius  of  the  cylinder  is  8.41  cm,  a  solution  plane  was 
specified  at  X  =  -9.19  and  the  spacing  on  either  side  of  this  plane  was  fixed  to  0.75  cm. 
Hence,  three  planes  are  specified  in  this  manner  starting  3.0  mm  from  the  leading  edge  of 
the  cyhnder  and  progressing  upstream  with  a  spacing  of  7.5  mm.  Also,  a  plane  was 
specified  at  the  unscoured  bed  level  (Z  =  0.0  cm)  with  a  spacing  of  5.0  mm.  For  the  case 
of  a  scour  hole,  an  additional  plane  was  specified  in  the  X-direction  at  the  leading  edge  of 
the  hole  in  order  to  better  capture  the  flow  separation.  For  all  cases,  the  mesh  extended 
from  -58.4  cm  to  -t40.0  cm  in  the  X-direction  and  -80.0  cm  to  0.0  cm  in  the  Y-direction. 
The  distance  in  the  Y-direction  dictates  that  the  spacing  between  piles  that  are  imposed 
on  the  solution  as  part  of  the  symmetry  boundary  condition  on  the  front  of  the  mesh  is  1.6 
m  from  center  to  center.  Thus,  the  piles  occupy  only  10%  of  the  cross-sectional  area 
between  the  piles  and  therefore  interference  between  piles  is  avoided.  The  distance  the 
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mesh  extended  in  the  Z-dkection  was  dependent  upon  the  depth  of  scour  hole  studied.  In 
each  direction  the  principle  followed  in  mesh  generation  was  that  spacing  between 
consecutive  planes  not  exceed  a  ratio  of  1.3.  This  criterion  is  suggested  for  solution 
stability.  For  each  of  the  simulations,  planes  were  added  or  subtracted  at  different  areas 
in  the  mesh  to  satisfy  this  rule.  The  largest  mesh  used  consisted  of  51,414  cells  and  the 
smallest  had  30,492  cells.  The  smallest  cell  generated  by  the  software  measured  0.75  x 
1.0x0.4  cm\ 


Front  and 
Back  BCs: 


Bottom  BC: 
Wall 


Figure  4.5  Summary  of  mesh  boundary  conditions. 
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For  the  initial  flat  bed  case,  the  only  obstacle  in  the  computational  mesh  required 
is  the  cylinder.  It  is  located  at  the  origin  in  the  mesh  and  extends  the  full  range  in  the  Z- 
direction.  It  has  a  diameter  of  16.8  cm  and  is  situated  on  the  back  face  of  the 
computationa:l  domain.  For  the  cases  involving  a  scour  hole,  a  much  more  complex 
obstacle  must  be  generated.  It  has  often  been  mentioned  in  the  Uterature  that  the  scour 
hole  resembles  an  inverted  right  cone  whose  axis  coincides  with  that  of  the  cylinder. 
Using  this  model,  the  scour  hole  is  constructed  by  first  creating  a  large  rectangular  block 
on  the  bottom  of  the  grid  and  then  excavating  out  a  hole  shaped  like  an  inverted  cone. 
Completion  of  the  model  is  achieved  by  next  adding  the  cylinder  at  the  center  of  the  hole. 
The  equation  for  the  cone  is  dependent  upon  both  the  dynamic  angle  of  repose  of  the 
sediment  and  the  desired  scour  depth.  This  technique  creates  a  sharp  interface  where  the 
cone  meets  the  cylinder.  Point  gauge  data  from  the  experiments  indicate  that  there  is  a 
narrow  but  noticeable  shelf  immediately  next  to  the  cylinder.  To  introduce  this  into  the 
model,  a  block  is  created  at  the  bottom  of  the  cone  to  truncate  the  hole  before  it  reaches 
the  cylinder.  For  simulations  of  the  experimental  work,  the  angle  of  the  cone,  width  of 
the  shelf  and  scour  depth  are  taken  dkectly  from  point  gauge  data.  Since  the  focus  of  this 
investigation  was  the  hydrodynamics  at  the  leading  edge  of  the  cylinder,  no  attempt  was 
made  to  model  the  characteristic  hump  that  is  found  in  the  wake  region  of  the  cylinder. 
The  combination  of  these  obstacles  creates  the  geometry  depicted  in  Figure  4.6. 


Cylinder 

Side  View  Top  View 

Figure  4.6  Scour  hole  obstacle  geometry. 
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4.2.2  Simulations  of  Fixed  Bed  Experiments 

The  four  physical  experiments  were  simulated  computationally  using  the 
techniques  described  in  the  previous  section.  Each  simulation  was  initiated  with  a  "hot- 
start"  and  run  for  20  seconds.  The  left  boundary  condition  was  determined  using  a  best 
fit  to  the  experimental  data  obtained  at  the  farthest  point  upstream  and  the  left  boundary 
of  the  mesh  was  set  at  that  point.  The  configuration  of  the  scour  hole,  especially  the 


depth  and  the  angle  of  the  cone,  was  determined  directly  from  the  point  gauge  data.  A 
summary  of  some  of  the  pertinent  simulation  parameters  is  found  in  Table  4. 1 . 


Tabl 

e  4,1  Physical  Experiment  Simulation  Parameters 

Scour  Hole  After 

Number  of 
Cells 

Smallest  Cell 
x,y,z  (cm') 

Largest  Cell 

x,y,2,  (cffl^) 

CPU  Time 
(hotirs) 

0  hours 

41,514 

0.75x1.0x0.5 

4.2x7.0x1.9 

32.67 

1  hour 

30,492 

0.75x1.0x0.5 

3.9x7.0x1.3 

15.75 

6  hours 

36,784 

0.75x1.0x0.5 

5.0x7.0x1.3 

15.46 

24  hours 

37,884 

0.75x1.0x0.4 

5.8x7.0x2.8 

13.48 

The  software  gives  a  number  of  graphs  as  standard  output.  They  include  graphs 
of  the  computational  grid,  the  fluid  fraction,  the  time  step  size,  the  number  of  pressure 
iterations,  the  volume  of  fluid  in  the  grid,  the  average  kinetic  energy  in  all  the  elements 
in  the  grid,  and  the  average  turbulent  kinetic  energy.  The  software  also  prepares  any 
graphs  requested  by  the  user.  For  the  simulations  in  this  section,  graphs  were  requested 
for  the  time  history  of  the  flow  properties  at  a  point  near  the  cylinder  at  an  elevation  equal 
to  the  unscoured  bed.  These  plots  indicate  whether  the  simulation  has  reached  an 
equilibrium  solution  at  the  end  of  the  simulation  time.  In  addition,  velocity  vector  plots 
with  pressure  contour  overlays  for  three  different  planes  were  requested.  Examples  of  the 
graphs  prepared  by  the  software  for  the  24  hour  scour  hole  are  contained  in  Figures  4.7 
(a)  through  (f). 
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Figure  4.7  Sample  output  from  FL0W3D®  for  the  simulation  of  the  scour  hole  after  24 
hours:  time  histories  of  (a)  X  velocity  at  the  base  of  the  cylinder,  (b)  estimated  mean 
kinetic  energy,  (c)  average  turbulent  kinetic  energy,  and  velocity  vector  plots  for  the 
planes:  (d)  Y  =  0  cm,  (e)  X  =  0  cm,  and  (f)  Y  =  -15.25  cm  (all  units  in  plots  are  cgs). 
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Figure  4.7  continued. 
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Figure  4.7  conlinucd. 
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From  Figures  4.7  (a),  (b)  and  (c),  it  is  apparent  that  the  model  has  indeed  reached 
an  equilibrium  state.  The  vector  plot  in  Figure  4.7  (d)  shows  the  clear  formation  of  the 
horseshoe  vortex.  In  order  to  compare  the  data  produced  from  the  simulations  to  the 
experimental  measurements,  it  is  necessary  to  interpolate  the  data  from  the  simulations  to 
the  points  were  the  measurements  were  made.  The  procedure  for  accompUshing  this  task 
involves  first  modifying  the  requested  output  file  to  direct  the  software  to  write  the  flow 
properties  (velocities,  turbulent  energy,  etc.)  to  an  ASCII  file.  This  file  is  then  read  and 
the  desired  velocities  at  the  points  where  measurements  are  made  are  found  through  a 
three  dimensional  interpolation  of  the  simulation  data.  The  program  compare.for  listed  in 
the  Appendix  was  used  for  this  feat.  Figures  4.8  through  4.14  contain  the  data  from  the 
simulations  displayed  in  the  same  manner  as  that  from  the  physical  experiments  in 
Chapter  3.  It  is  important  to  note  that  the  bottom  boundary  in  the  graphs  is  from  the  point 
gauge  data  for  the  physical  experiments  and  is  not  the  bottom  boundary  used  in  the 
simulations.  For  the  following  plots  the  velocity  vector  length  is  proportional  to  its 
magnitude  and  it  is  scaled  by  the  magnitude  of  the  largest  vector  found  in  the  plane.  This 
magnitude  is  listed  in  all  plots.  Also,  the  contours  show  the  magnitude  of  the  velocity 
normal  to  the  plane  with  a  positive  magnitude  indicating  flow  into  the  page.  The 
coordinate  system  is  the  same  as  described  in  Chapter  3.  Figures  4.8,  to  4.1 1  show  the 
velocities  for  the  scour  holes  at  times  t  =  0,  1,6  and  24  hours  respectively.  Figure  4.12 
shows  the  bottom  velocities  at  the  same  elevation  as  those  in  Chapter  3.  Care  should  be 
taken  in  interpreting  the  results  shown  in  this  figure  as  not  all  the  measurements  were 
taken  at  the  same  elevation  off  the  bed.  Figures  4.13  and  4.14  show  the  velocities  in  the 
planes  Z  =  +1  cm  and  Z  =  +5  cm  where  Z  =  0  cm  is  the  level  of  the  unscoured  bed. 
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Figure  4.8  Mean  Velocities  for  the  initial  flat  bed  case  in  the  planes:  (a)  0°,  (b)  45°,  (c) 
90°,  (d)  Y=1.8  cm,  (e)  Y=4.3  cm,  (f)  Y=7.3  cm,  (g)  X=-b,  and  (h)  X=-2b. 
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Figure  4.8  continued. 
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Figure  4.8  continued. 
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Figure  4.8  continued. 
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T  =  1  Hours,  0  Degrees 
Simulation  Output 


Figure  4.9  Mean  Velocities  for  the  t  =  1  hour  scour  hole  in  the  planes: 
(a)  0°,  (b)  30°,  (c)  60°,  (d)  90°,  (e)  X=-b,  and  (f)  X=-2b. 
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Figure  4.9  continued. 
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T  =  1  Hours,  X  =  -b 
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Figure  4.9  continued. 
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T  =  6  Hours,  0  Degrees 
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Figure  4.10  Mean  Velocities  for  the  t  =  6  hours  scour  hole  in  the  planes: 
(a)  0°,  (b)  30°,  (c)  60°,  (d)  90°,  (e)  X=-b,  and  (t)  X=-2b. 
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Figure  4.10  continued. 
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Figure  4.10  continued. 
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T  =  24  Hours,  0  Degrees 
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Figure  4.1 1  Mean  Velocities  for  the  t  =  24  hour  scour  hole  in  the  planes: 
(a)  0°,  (b)  30°,  (c)  60°,  (d)  90°,  (e)  X=-b,  and  (f)  X=-2b. 


Figure  4. 1 1  continued. 
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Figure  4.1 1  continued. 


Figure  4.12  Mean  Velocities  at  the  bed  at  times:  (a)  t  =  0  hours, 
(b)  t  =  1  hour,  (c)  t  =  6  hours,  and  (d)  t  =  24  hours. 
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Figure  4.12  continued. 


Figure  4.13  Mean  Velocities  for  the  plane  Z  =  +1  cm  at  times:  (a)  t  =  0  hours, 
(b)  t  =  1  hour,  (c)  t  =  6  hours,  and  (d)  t  =  24  hours. 


Figure  4.13  continued. 


Figure  4. 14  Mean  Velocities  for  the  plane  Z  =  +5  cm  at  times:  (a)  t  =  0  hours, 
(b)  t  =  1  hour,  (c)  t  =  6  hours,  and  (d)  t  =  24  hours. 


Figure  4.14  continued. 
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A  thorough  discussion  of  the  results  of  the  simulations  and  their  comparison  with 
the  experimental  measurements  is  presented  in  Chapter  5.  The  following  chapter 
contains  error  analyses  of  both  the  magnitude  and  direction  of  the  velocities  predicted  by 
the  software.  The  preceding  figures  confirm  a  number  of  the  observations  made  about  the 
experimental  data.  The  velocity  vectors  for  the  flat  bed  case  (Figure  4.8)  again  depict  the 
lack  of  a  horseshoe  vortex.  The  main  hydrodynamic  mechanisms  in  evidence  are  the 
downflow  at  the  leading  edge  and  the  acceleration  of  the  flow  around  the  sides  of  the  pile. 
The  acceleration  to  the  sides  is  especially  evident  in  the  series  of  graphs  where  Y  =  0  cm, 
1.8  era,  4.3  cm  and  7.3  cm.  The  contours  of  the  normal  velocity  show  a  steady  increase 
in  magnitude  as  one  moves  away  from  the  stagnation  plane. 

Figure  4.9  shows  the  scour  hole  after  1  hour  of  scouring.  Although  there  appears 
to  be  more  vorticity  at  the  very  bottom  of  the  hole,  there  is  as  yet  no  definable  horseshoe 
vortex.  There  is  also  no  separation  at  the  leading  edge  of  the  hole.  This  is  well  illustrated 
by  graphs  showing  the  0°  and  30°  planes.  The  60°  and  90°  planes,  however,  start  to  show 
some  of  the  counter-clockwise  rotation  that  is  associated  with  this  vortex. 

Figure  4.10  contains  graphs  of  the  mean  velocities  after  6  hours  of  scouring. 
Here,  the  separation  of  the  flow  at  the  leading  edge  of  the  pile  is  readily  apparent.  The 
horseshoe  vortex  is  quite  well  defined  and  occupies  almost  the  entire  hole  cross-section  in 
the  stagnation  plane.  As  the  vortex  is  swept  around  the  sides  of  the  pile  the  vortex  seems 
to  narrow  and  is  swept  up  the  sides  of  the  hole.  There  is  also  an  intensification  of  the 
normal  velocity  towards  the  center  of  the  vortex. 

Figure  4.1 1  shows  the  flow  field  around  the  scour  hole  as  the  equilibrium  value  is 
neared.  The  separation  of  the  flow  at  the  leading  edge  of  the  hole  is  much  more 
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pronounced.  The  horseshoe  vortex  occupies  the  entire  hole  in  the  stagnation  plane  and 
narrows  rapidly  as  it  is  swept  around  the  edges.  The  movement  of  the  vortex  up  the  sides 
of  the  hole  as  it  is  swept  around  is  also  more  apparent  in  this  case  than  for  that  of  the  6 
hour  hole. 

The  intentional  lack  of  care  in  modeling  the  wake  region  makes  the  results  there 
suspect.  However,  Figures  4.8  through  4.1 1  show  the  same  behavior  as  the  experimental 
measurements.  There  is  a  low  velocity  region  immediately  behind  the  cylinder  and  a  well 
defined  upwelling.  The  series  of  figures  for  X  =  -b  depict  a  jump  in  the  magnitude  of  the 
velocity  in  the  plane  from  the  flat  bed  case  to  the  t  =  1  hour  case.  After  this  maximum, 
the  magnitude  of  the  largest  vector  decreases  for  the  t  =  6  hour  and  t  =  24  hour  cases. 
This  is  exactly  the  same  behavior  found  in  the  physical  experiments  though  less 
pronounced. 

The  bottom  velocities  are  shown  in  Figure  4.12.  There  is  a  very  clear  difference 
in  the  flow  patterns  as  the  scour  hole  deepens.  Initially,  the  pattern  is  dominated  by  the 
acceleration  of  flow  around  the  sides  of  the  pile  as  is  seen  in  graphs  (a)  and  (b).  As  time 
elapses,  a  horseshoe  vortex  appears  and  it  becomes  the  dominant  pattern  at  the  leading 
edge  as  evidenced  by  the  reversal  of  flow  direction  in  the  stagnation  plane  in  graphs  (c) 
and  (d).  Figures  4.13  and  4.14  show  the  expected  two  dimensional  acceleration  of  the 
flow  around  the  sides  of  the  cylinder.  The  growth  of  the  magnitudes  of  the  contours  show 
the  flow  into  and  out  of  the  scour  hole  as  it  deepens. 
4.2.3  Simulation  of  the  Scour  Hole  as  a  Function  of  Time 

In  order  to  study  the  hydrodynamics  of  the  scour  hole  as  a  function  of  time,  a 
more  comprehensive  data  set  covering  the  evolution  of  the  hole  must  be  obtained.  To  this 
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end,  a  series  of  simulations  was  performed  with  the  same  scour  situation  as  the  physical 
experiments.  Since  the  scour  depth  is  a  function  of  time,  the  approach  taken  was  to 
perform  a  number  of  simulations  at  different  depths  as  the  scour  hole  grows  and  analyze 
the  evolution  of  the  hydrodynamics  associated  with  the  different  depths.  The  series 
modeled  the  scour  hole  at  depths  occurring  at  0.5  cm  and  1.0  cm  intervals  from  the  flat 
bed  case  to  the  24  hour  case.  A  total  of  21  simulations  was  run.  The  hole  was 
constructed  in  the  same  manner  as  described  in  section  4.2.1.  The  angle  of  repose  of  the 
sediment  was  determined  using  the  average  angle  calculated  from  the  point  gauge  data  for 
the  stagnation  plane  for  the  four  cases.  A  value  of  (j)  =  31.2°  was  used.  The  left,  user 
specified  boundary  condition  was  found  using  a  multiple  least  squares  fit  to  the  data 
occurring  the  farthest  upstream  from  the  physical  experiments  for  all  four  cases.  The 
velocity  profile  used  is  illustrated  in  Figure  4.15. 

The  simulation  times  for  the  series  was  20  seconds  and  the  simulations  were  all 
initialed  with  a  "hot-start."  The  turbulence  closure  scheme  used  was  the  RNG  turbulence 
model.  In  addition  to  the  28  cases,  a  threshold  condition  was  run.  This  simulation 
corresponds  to  the  depth  averaged  velocity  ratio  where  scour  initiates,  U/Uc  =  0.45.  The 
left  boundary  condition  for  this  case  is  simply  the  velocity  profile  for  the  series  multiplied 
by  U/Uc  for  the  threshold  condition  and  divided  by  UAJc  for  the  physical  experiments. 
Presentation  of  the  results  for  these  simulations  and  an  analysis  of  the  evolution  of  scour 
is  contained  in  Chapter  5. 
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Figure  4.15  Best  fit  velocity  profile  to  the  experimental  data. 


CHAPTER  5 
RESULTS  AND  DISCUSSION 


This  chapter  integrates  the  physical  experiments  and  the  numerical  simulations 
performed  as  part  of  this  investigation.  The  chapter  is  divided  into  two  sections.  The 
first  section  describes  the  ability  of  the  computational  model  to  reproduce  the  velocity 
measurements  made  during  the  physical  experiments.  The  second  section  discusses  a 
new  procedure  for  using  the  FL0W-3D®  software  to  predict  the  scour  depth  as  a  function 
of  time. 

5.1  Numerical  Simulation  Accuracy 
The  computational  model  is  judged  for  its  accuracy  in  two  ways:  its  ability  to 
predict  the  velocity  magnitude  and  its  abiHty  to  determine  the  velocity  direction.  As 
described  in  the  previous  chapter,  in  order  to  compare  the  output  of  FL0W-3D®  with  the 
measurements,  a  three-dimensional  interpolation  of  the  model  was  performed  to 
determine  the  velocity  at  the  same  point  where  the  measurements  were  taken.  This  task 
was  accomplished  by  the  program  compare.for  found  in  the  Appendix.  The  error  in  the 
magnitude  of  the  velocity  is  represented  by  the  percentage  error  given  in  equation  (5.1). 
In  equation  (5.1),  the  subscript  e  indicates  the  value  from  the  experiment  and  the 
subscript  c  indicates  the  result  from  the  computational  simulation. 
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%  Error  in  Velocity  Magnitude  = 


2,2,  2^ 


xlOO     (5. 1) 


The  error  in  velocity  direction  is  represented  by  the  three-dimensional  angle  between  the 
computed  and  measured  velocity  vectors  (in  degrees).  This  angle  is  found  by  taking  the 
inverse  cosine  of  the  normalized  dot  product  between  the  measured  and  computed 
vectors.  The  results  of  this  analysis  can  be  summarized  by  the  histograms  in  Figures  5.1 
through  5.5. 
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78%  of  measurements  have  less 
than  30%  difference  in  magnitude. 
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81%  of  measurements  have  less 
than  30  degrees  difference  in  angle. 
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Figure  5.  1  Histograms  of  the  total  error  from  all  the  experiments  (3906  measurements). 
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Figure  5.  2  Histograms  of  error  for  the  T  =  0  hours  experiment  (unscoured  bed). 

The  histogram  in  Figure  5.1  shows  the  total  error  for  all  the  3906  measurements 
performed.  Each  of  the  bars  for  the  magnitude  percentage  error  in  Figures  5.1  through 
5.5  represents  the  number  of  measurements  having  the  indicated  error  at  10%  increments. 
Each  of  the  bars  for  the  error  in  angle  represents  the  number  of  measurements  having  the 
indicated  error  at  10  degree  increments.  Overall,  the  model  performed  reasonably  well. 
It  predicted  the  magnitude  of  the  velocity  within  30%  error  for  78%  of  the  measurements 
and  the  direction  within  30  degrees  for  81%  of  the  measurements. 
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T  =  1  Hour  Experiment 
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69%  of  measurements  have  less 
than  30%  difference  in  magnitude. 


I  fl — LLlJJ  1 1  I  1 1   I  n~\ I — ^ m        ■   I  I  I— I— > 

20       40        60       80       100      120      140      160      180  200 
Magnitude  Error  in  Percent  Difference 

T  I  1  1  1  1  1  1  


89%  of  measurements  have  less 
than  30  degrees  difference  in  angle. 
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Figure  5.  3  Histograms  of  error  for  the  T  =  1  hour  experiment. 

For  each  of  the  figures,  the  200%  bin  for  the  magnitude  error  appears 
unacceptably  large.  Since  this  is  the  cut-off  for  the  histogram,  this  bin  shows  the  number 
of  experiments  having  a  percentage  difference  greater  than  190%.  The  reason  that  the 
error  in  magnitude  can  be  this  great  is  due  to  the  fact  that  at  the  boundaries  (bed  and 
cylinder)  the  measured  velocities  are  quite  low.  Therefore,  since  the  percentage  error  is 
normalized  by  this  value,  very  large  percentages  can  occur.  For  example,  if  the  measured 
value  is  0.5  cm/s  and  the  computed  value  is  2  cm/s  a  300%  error  results. 
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78%  of  measurements  have  less 
than  30%  difference  in  magnitude. 
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78%  of  measurements  have  less 
than  30  degrees  difference  in  angle. 
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Figure  5.  4  Histograms  of  error  for  the  T  =  6  hour  experiment. 
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Individually,  the  model  for  the  24  hour  scour  hole  performed  the  best  in  the 
prediction  of  the  velocity  magnitude  and  the  model  for  the  one  hour  scour  hole  performed 
the  best  in  predicting  the  direction  of  the  velocity.  Conversely,  the  model  for  the  one 
hour  scour  hole  performed  the  worst  in  the  prediction  of  the  velocity  magnitude  and  the 
worst  predictor  of  the  velocity  direction  was  the  model  for  the  six  hour  scour  hole 
followed  closely  by  the  model  for  the  24  hour  scour  hole.  This  indicates  an  inverse 
relationship  in  the  models  ability  to  predict  both  the  magnitude  and  the  direction 
accurately. 
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87%  of  measurements  have  less 
than  30%  difference  in  magnitude. 
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79%  of  measurements  have  less 
than  30  degrees  difference  in  angle. 
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Figure  5.  5  Histograms  of  error  for  the  T  =  24  hours  experiment. 

Contour  plots  of  the  error  in  velocity  magnitude  and  difference  in  direction  were 
made  for  each  of  the  vertical  planes  presented  in  chapters  3  and  4.  In  addition,  plots  of 
the  difference  vectors  were  constructed.  The  vectors  in  these  plots  show  the  experimental 
velocity  vector  minus  the  computational  velocity  vector  in  the  plane  depicted  and  the 
contours  show  the  difference  in  the  velocity  magnitude  (in  cm/s)  in  the  direction 
perpendicular  to  the  plane.  Due  to  the  large  number  of  plots  created  (78)  only  a  few  wiU 
be  presented  in  this  chapter  to  illustrate  the  model's  accuracy.  All  of  the  plots  showing 
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the  error  and  difference  vectors  are  contained  in  the  reference  Gosselin  and  Sheppard 
(1997). 

The  model  performed  best  in  the  front  quadrant  of  the  cylinder  and  the  worst  in 
the  wake  region.  This  result  is  expected  since  care  was  taken  to  model  the  scour  hole  in 
the  front  and  less  care  was  taken  in  modeling  the  geometry  of  the  bed  in  the  lee  of  the 
structure.  In  the  stagnation  plane,  the  model  did  well  in  predicting  magnitude  but 
relatively  poorly  in  prediction  of  the  direction.  Figure  5.6  and  5.7  show  the  scour  hole 
after  6  hours.  While  the  majority  of  the  error  is  less  than  30%  in  the  velocity  magnitude, 
the  difference  in  angle  reaches  values  upwards  of  75  degrees. 
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Figure  5.6  Results  for  the  6  hour  scour  hole.  The  difference  in  velocity  vectors  in 
the  stagnation  plane  (0  =  0°)  are  depicted  by  arrows.  The  contours  are  differences 
in  velocity  (in  cm/s)  in  the  direction  normal  to  the  stagnation  plane. 


Figure  5.7  Contours  of  error  in  (a)  velocity  magnitude  and  (b)  direction  in 
the  stagnation  plane  (6  =  0°)  for  the  6  hour  scour  hole. 
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Most  of  the  error  in  velocity  direction  occurs  within  the  scour  hole.  Again,  this 
result  is  expected.  Since  the  condition  for  the  stagnation  plane  in  the  computer  model  is  a 
symmetry  plane,  there  can  be  no  velocity  component  normal  to  this  plane  (see  Figure 
4.10(a)).  In  the  physical  experiment,  the  existence  of  a  true  vertical  stagnation  plane  is 
unlikely.  Also,  the  Acoustic  Doppler  Velocimeter  averages  the  velocity  vector  over  a 
9ram  long  finite  control  volume  as  opposed  to  an  infinitesimal  point.  These  combine  to 
make  measurement  of  a  finite  flow  normal  to  the  stagnation  plane  highly  likely.  Since  the 
measured  velocities  are  low  near  the  cylinder,  any  flow  normal  to  the  plane  results  in  a 
significant  change  in  du-ection  of  the  velocity  vector.  Hence,  error  in  vector  du-ection  is 
large  in  the  stagnation  plane  near  the  cylinder  and  in  the  scour  hole  for  all  cases. 


T  -  24  Hours,  90  Degrees 
Difference  Vectors  and  Contours 


Radial  Position  (cm) 

Figure  5.8  Results  for  the  24  hour  scour  hole.  The  difference  in  velocity  vectors 
in  the  90°  plane  (0  =  90°)  are  depicted  by  arrows.  The  contours  are  differences  in 
velocity  (in  cm/s)  in  the  direction  normal  to  the  stagnation  plane. 


Figure  5.9  Contours  of  error  in  the  (a)  velocity  magnitude  and  (b) 
direction  in  the  90°  plane  (0  =  90°)  for  the  24  hour  scour  hole. 
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The  model  performed  particularly  well  in  the  30°,  60°  and  90°  planes.  Figure  5.8 
and  5.9  show  the  difference  vectors  and  errors  for  the  plane  located  90°  from  the 
stagnation  plane  (at  X  =  0.0  cm)  for  the  24  hour  scour  hole.  These  figures  indicate  that 
the  majority  of  the  vectors  are  predicted  within  30%  of  the  velocity  magnitude  and  within 
20°  of  the  velocity  direction.  These  figures  are  typical  of  all  the  planes  other  than  the 
stagnation  or  the  wake.  Not  coincidentally,  these  planes  also  contain  the  areas  where  the 
measured  velocities  are  the  largest.  The  larger  velocities  lead  to  more  forgiveness  in  the 
estimation  of  the  error  in  velocity  magnitude. 


T  =  1  Hours,  X  a-bi 
Difference  Vectors  and  Contours 


j  largest  vector  =  12.6cfn/s 

~^-5  0  5  10  15  20  25  30 

Y  Position  (cm) 

Figure  5.10  Results  for  the  one  hour  scour  hole.  The  difference  in  velocity  vectors 
in  the  wake  region  (X  =  -b)  are  depicted  by  arrows.  The  contours  are  differences 
in  velocity  (in  cm/s)  in  the  direction  normal  to  the  stagnation  plane. 
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Figure  5.11  Contours  of  error  in  the  (a)  velocity  magnitude  and 
(b)  direction  in  the  wake  region  (X  =  -b)  for  the  one  hour  scour  hole. 
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Finally,  the  model  performed  poorest  in  the  planes  measured  behind  the  cylinder 
at  the  distances  X  =  -b  and  X  =  -2b.  Figure  5.10  and  5.11  show  the  error  and  difference 
vectors  for  the  plane  X  =  -b  for  the  one  hour  scour  hole.  As  mentioned  previously,  this 
result  is  expected  since  the  wake  region  was  not  modeled  accurately.  The  accuracy  of  the 
model  decreases  sharply  directly  behind  the  cyhnder  (Y  >  8.4cm).  The  error  in  modeling 
the  wake  region  improves  both  with  increased  scour  depth  and  with  distance  from  the 
cylinder. 

In  conclusion,  the  models  were  deemed  capable  of  predicting  the  velocity 
magnitude  and  direction  within  a  satisfactory  margin  of  error.  It  was  therefore  possible  to 
continue  with  the  part  of  the  investigation  involving  the  prediction  of  scour  as  a  function 
of  time. 

5.2  Prediction  of  Scour  Depth  as  a  Function  of  Time 
The  crux  of  this  investigation  is  to  develop  a  method  for  examining  the  depth  of 
scour  as  a  function  of  time.  To  this  end,  a  series  of  simulations  was  run  using  average 
values  from  the  physical  experiments  as  input.  The  simulations  modeled  the  physical 
experiment  at  scour  depths  ranging  from  0.0cm  to  13.5cm  at  intervals  equal  to  0.5cm  and 
1.0  cm.  The  conditions  for  the  simulations  are  summarized  in  Table  5.1.  A  total  of  21 
simulations  was  run  including  a  simulation  over  a  flat  bed  initial  condition  at  the 
threshold  velocity  ratio  (U/Uc  -  0.45).  From  the  simulations  the  velocity  and  tufbulent 
kinetic  energy  data  were  gathered  in  ASCII  format. 


Table  5.1  Summary  of  simulations  for  determining  scour  depth  as  a  function  of  time. 

Depth  Averaged  Velocity,  U 

25.74  cm/s 

Water  Depth,  yo 

34.86  cm 

Cylinder  Diameter,  b 

16.83  cm 

Dynamic  Angle  of  Repose,  (j) 

31.22° 
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Since  there  is  a  wealth  of  hydrodynamic  data  available  from  the  numerical 
simulations,  a  reasonable  approach  to  predicting  the  scour  depth  time  history  would 
involve  relating  the  hydrodynamics  inside  the  scour  hole  to  the  rate  of  scour.  This 
approach  makes  good  physical  sense  since  it  is  the  hydrodynamics  which  drive  the 
sediment  transport.  Examining  the  sediment  transport  of  the  scour  event,  for  clear  water 
scour,  the  important  parameter  is  the  outflux  of  sediment  from  the  scour  hole,  Qo.  Qo 
varies  significantly  with  the  scour  depth.  Initially  it  is  relatively  large,  scouring  the  hole 
out  rapidly.  It  then  tapers  off  and  approaches  zero  as  the  hole  approaches  its  equilibrium 
depth.  Employing  the  same  assumption  used  in  the  numerical  simulations,  the  hole  is 
modeled  as  an  inverted  right  cone  whose  axis  coincides  with  that  of  the  cylinder.  The 
cone  is  truncated  such  that  a  ledge  of  a  known  width  is  formed  at  the  specific  scour  depth. 
For  the  conditions  of  the  experiments  performed  at  the  University  of  Florida,  it  has  been 
observed  that  the  width  of  the  ledge  formed  in  the  scour  hole  is  on  the  order  of  1  cm. 
This  width  has  proved  to  be  relatively  insensitive  to  both  sediment  size  and  pile  diameter 
for  the  range  of  values  tested. 

The  sediment  outflux  is  related  to  the  rate  of  change  of  the  volume  of  the  cone. 
Including  the  effects  of  the  volume  concentration  of  the  sediment,  Qo  can  be  expressed  as 
a  function  of  the  scour  depth  and  the  rate  of  change  of  the  scour  depth  as  shown  in 
equation  (5.2).  In  this  equation,  V  is  the  volume  and  ds  is  the  depth  of  the  scour  hole,  rc 
is  the  radius  of  the  cyhnder,  and  wl  is  the  width  of  the  ledge.  In  equation  (5.2),  c  is  the 
volume  concentration  of  the  sediment  (-  volume  of  the  grains  /  bulk  volume). 
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dV  d  f  Tidl  ,  <(r,+wj 
—  —  f — <  
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dt        dt[  Stance])  taiKj) 


(5.2) 


Therefore,  if  Qo  is  a  known  function  of  scour  depth,  then  the  time  history  of  scour  depth 
can  be  calculated  by  solving  a  nonlinear  first  order  differential  equation  of  the  form 


=  f(dj. 


The  problem  hinges  on  finding  a  method  for  determining  Qo  from  the 
hydrodynamic  properties.  The  most  likely  candidate  for  this  should  be  the  shear  stress  at 
the  bed.  In  general,  however,  an  accurate  calculation  of  bed  shear  stress  is  difficult.  This 
is  particularly  true  when  using  FL0W-3D®.  Since  the  grid  is  orthogonal  rather  than 
body-fitted,  the  points  where  the  velocities  are  calculated  are  found  at  varying  distances 
from  the  bottom  in  the  scour  hole.  Also,  while  the  calculation  of  the  shear  stress  under 
open  channel  flow  over  a  flat  bed  is  well  estabhshed,  the  methodology  for  calculation  of 
the  shear  stresses  in  the  scour  hole,  where  both  the  geometry  is  more  complex  and 
swirling  /  impinging  flows  are  found,  is  less  obvious.  For  these  reasons,  an  alternate  flow 
property  was  sought  as  an  independent  variable  for  Qo. 

In  order  to  find  an  appropriate  property,  it  is  helpful  to  examine  the 
hydrodynamics  required  to  erode  a  single  sediment  particle  from  the  scour  hole.  The 
particle  must  first  be  subjected  to  a  sufficient  lift  and  drag  force  to  dislodge  it  from  its 
resting  point  at  the  bed.  Since  it  is  the  flow  that  is  doing  work  on  the  bed  in  the  form  of 
sediment  transport,  an  examination  of  the  energy  associated  with  the  flow  near  to  the  bed 
seems  in  order.  At  the  initial  stages  of  scour,  a  particle  that  is  given  a  kinetic  energy 
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above  a  certain  threshold  value  can  be  removed  from  the  proximity  of  the  piling  without  a 
significant  vertical  displacement;  i.e.,  it  can  be  transported  away  in  the  form  of  bed  load. 
Once  the  scour  hole  has  formed,  however,  a  particle  must  also  overcome  a  potential 
energy  related  to  the  depth  of  scour.  Using  this  argument,  total  flux  of  sediment  out  of 
the  hole,  Qo,  should  be  related  to  both  the  excess  kinetic  energy  at  some  (short)  fixed 
distance  above  the  surface  of  the  hole  and  the  potential  energy  associated  with  removing 
sediment  from  the  hole.  The  use  of  the  excess  kinetic  energy  (i.e.  the  kinetic  energy 
minus  the  threshold  kinetic  energy)  incorporates  the  fact  that  a  certain  level  of  energy 
must  be  exceeded  in  order  to  remove  the  sediment. 

To  account  for  the  fact  that  sediment  may  be  dislodged  at  any  point  on  the  surface 
of  the  hole,  the  integral  of  the  excess  kinetic  energy  over  the  entire  surface  is  examined. 
Once  a  particle  is  dislodged,  it  enters  the  flow  field  in  the  hole.  One  of  the  dominant 
features  of  the  flow  field  towards  the  front  of  the  cylinder  is  the  downwards  flow 
resulting  from  the  stagnation  pressure  gradient.  Thus  the  particle  in  the  flow  will  most 
likely  be  transported  near  the  bottom  of  the  hole  before  it  is  removed.  Therefore,  the  use 
of  the  deepest  point  of  scour  (depth  of  scour)  appears  to  be  a  good  choice  for  a 
representative  value  of  the  potential  energy  that  must  be  overcome. 

The  choices  of  these  two  flow  properties  must  be  reconciled  with  what  is  known 
about  the  equilibrium  scour  depth  since  that  is  the  state  that  is  approached  through  the 
temporal  variation  of  local  scour.  Many  researchers  (e.g.  Sheppard  (1996)  and  Chiew 
(1984))  have  documented  that  under  conditions  of  clear  water  scour,  the  equilibrium 
depth  is  a  strong  function  of  the  pile  diameter  (b),  upstream  water  depth  (yo),  sediment 
diameter  (dso),  and  depth  averaged  velocity  ratio  (U/Uc).  These  parameters  must  also 
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influence  the  properties  chosen  for  this  analysis.  The  influence  of  cylinder  diameter 
shows  up  in  the  integral  of  the  excess  kinetic  energy.  The  larger  the  diameter,  the  greater 
the  surface  area  of  the  hole  for  a  given  scour  depth.  A  larger  surface  means  a  greater  area 
over  which  the  kinetic  energy  may  be  above  the  threshold  value  and  hence  a  larger  Qo. 
The  upstream  water  depth  affects  the  kinetic  energy  at  the  bed  surface  only  when  it  is 
sufficiently  shallow.  It  is  only  when  the  presence  of  the  free  surface  is  "felt"  at  the  bed 
level  that  this  parameter  becomes  significant  in  the  determination  of  kinetic  energy.  This 
is  consistent  with  what  is  known  about  the  aspect  ratio  (yo/b,  see  Chapter  2). 

The  sediment  diameter  affects  the  kinetic  energy  in  a  number  of  ways.  First,  the 
sediment  diameter  is  directly  related  to  the  roughness  of  the  bed.  It  is  the  roughness, 
along  with  the  depth  averaged  velocity,  that  defines  the  shape  of  the  upstream  velocity 
profile.  This  profile  has  a  major  impact  on  the  secondary  flows  (the  horseshoe  vortex  and 
the  downflow  at  the  leading  edge)  inside  the  hole.  These  parameters  will  obviously  affect 
the  kinetic  energy  near  the  bed.  Next,  the  dynamic  angle  of  repose  of  the  sediment  is  a 
function  of  sediment  size.  The  angle  of  repose  will  affect  the  shape  of  the  hole  and  hence 
both  the  hydrodynamics  and  the  limits  of  the  integral.  Finally,  the  sediment  diameter  and 
mass  density  determine  the  critical  shear  stress  and  thus  the  threshold  kinetic  energy 
required.  The  final  property,  depth  averaged  velocity  ratio,  has  an  obvious  influence  on 
the  kinetic  energy  since  the  energy  is  directly  related  to  the  square  of  the  velocity. 

The  threshold  kinetic  energy  can  be  calculated  by  analyzing  the  threshold  scour 
condition.  It  is  generally  agreed  upon  in  the  literature  that  local  scour  around  a  circular 
cylinder  initiates  at  a  depth  averaged  velocity  ratio  (UAJc)  approximately  equal  to  0.45. 
This  is  the  threshold  scour  condition.  The  maximum  kinetic  energy  found  in  this 
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situation  should  give  the  threshold  value  for  kinetic  energy.  In  other  words,  it  is  assumed 
that  if  a  value  of  kinetic  energy  (KE)  found  inside  the  scour  hole  is  greater  than  the 
threshold  kinetic  energy  (KEt)  then  scour  is  still  occurring.  With  further  analysis,  it  is 
reasonable  to  assume  that  the  threshold  kinetic  energy  at  some  (short)  fixed  distance  from 
the  bed  could  be  determined  in  the  same  manner  as  the  critical  depth  averaged  velocity. 

The  total  kinetic  energy  per  unit  mass  used  for  the  analysis  is  given  by  equation 
5.3,  where  k  is  the  turbulent  kinetic  energy. 

Total  Kinetic  Energy  =  KE  =  ^(u' +  v' +  w')  +  k  (5.3) 

The  outputs  from  the  simulations  were  reduced  to  yield  the  KE  at  a  perpendicular 
distance  of  5.0  mm  from  the  bed.  The  choice  of  5.0  mm  reflects  the  limit  of  reliable 
measurements  in  the  physical  model  and  thus  rehable  comparisons  between  the 
measurements  in  the  physical  experiments  and  the  numerical  simulations.  The  KE  was 
determined  over  a  dense  orthogonal  mesh  with  a  5.0  mm  spacing  over  a  "surface"  5.0  mm 
above  the  scoured  bed.  The  velocities  for  the  calculation  were  found  using  a  three 
dimensional  interpolation  scheme.  This  procedure  is  contained  in  the  program 
"energy.for"  located  in  the  Appendix.  Examples  of  the  results  given  by  the  model  are 
illustrated  in  the  contour  plots  in  Figures  5.12  through  5.15.  These  figures  depict  the  hole 
at  depths  approximately  equal  to  those  in  tlie  physical  experiments. 

Figure  5. 1 6  shows  the  KE  over  the  bed  for  the  threshold  condition  case.  This  case 
was  run  using  a  similariy  shaped  velocity  profile  (see  Figure  4.15)  as  the  rest  of  the 
simulations  but  with  a  depth  averaged  velocity  ratio  equal  to  0.45.  For  the  threshold 
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condition,  the  maximum  KEt  found  has  a  value  of  125  cm^/s^  and  is  located  at  a  radial 
distance  of  1 .7  cm  from  the  cylinder  at  an  angle  of  70°  from  the  stagnation  plane. 


Total  Kinetic  Energy  Contours  (d  =  0.0  cm) 
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Figure  5.12  Total  kinetic  energy  contours  for  the  flat  bed  (ds  =  0.0  cm)  case. 


Total  Kinetic  Energy  Contours  (d  =  6.0  cm) 
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Figure  5.13  Total  kinetic  energy  contours  for  the  ds  =  6.0  cm  case. 


178 


Total  Kinetic  Energy  Contours  (d^  =  11.0cm) 
KE  =  (l/2)(u^  +      +  w^)  + k 
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Figure  5.14  Total  kinetic  energy  contours  for  the  ds  =  1 1.0  cm  case. 


Total  Kinetic  Energy  Contours  (d^  =  13.5  cm) 
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Figure  5.15  Total  kinetic  energy  contours  for  the  ds  =  13.5  cm  case. 
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Total  Kinetic  Eiiergy  Contours  (Threshold  Cond.,  d  =  0.0  cm) 
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Figure  5.16  Total  kinetic  energy  contours  for  the  threshold  condition  case. 

The  proposed  hypothesis  is  that  the  volumetric  flux  of  sediment  out  of  the  scour 
hole  is  proportional  to  the  excess  KE  over  the  entire  area  of  the  hole  and  the  potential 
energy  required  to  remove  sediment  from  the  bottom  of  the  hole.  Symbolically,  this 
concept  is  illustrated  in  equation  5.4. 

Q„  =       jj(KE-KE,)da  -  C,(gdJ.  (5.4) 

hole 

In  this  equation,  Ci  and  Cz  are  dimensional  coefficients  of  proportionality.  The 
coefficients  themselves  are  functions  of  the  flow  or  sediment  parameters.  The  equation 
shows  the  competing  processes  of  the  erosive  and  resistive  actions  of  the  hydrodynamics 
and  the  geometry  of  the  hole.  The  first  term  on  the  right  hand  side  represents  the  ability 
of  the  hydrodynamics  to  remove  sediment.  It  incorporates  into  it  all  the  factors  affecting 
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the  velocity  of  the  fluid  inside  the  hole.  As  a  result,  it  is  suspected  that  the  coefficient  Ci 
is  most  likely  a  function  of  such  parameters  as  U/Uc,  yo/b,  and  dso/b.  The  second  term  on 
the  right  hand  side  represents  the  barrier  to  be  overcome  in  order  to  remove  sediment 
from  the  hole.  This  term  is  a  function  of  fee  geometry  of  the  hole  and  effective  density  of 
the  sediment.  Removal  of  sediment  not  only  requires  a  vertical  displacement,  but  also  a 
translation  horizontally,  thus  C2  may  also  be  a  function  of  the  fall  velocity  of  the  sediment 
particle.  Note  that  the  effect  of  the  sediment  porosity  can  be  absorbed  into  both  of  the 
coefficients.  The  form  of  this  equation  is  also  attractive  because  it  clearly  defines  the 
equilibrium  condition.  Sediment  outflux  becomes  zero  when  the  energy  eroding 
sediment  is  equal  to  the  energy  necessary  to  remove  a  particle  from  the  hole. 

To  test  this  hypothesis,  the  equation  derived  for  Qo  based  on  the  rate  of  change  of 
the  volume  of  the  scour  hole  (equation  (5.2),  c  =  0.63)  was  compared  with  the  equation 
just  proposed  (equation  (5.4)).  Both  were  plotted  versus  scour  depth  and  the  coefficients 
Ci  and  C2  were  found  using  a  least  squares  fit  TTie  rate  of  scour  depth  term  in  equation 
'd(ds)^ 


(5.2), 


d.  =A 


dt 


,  was  determined  by  fitting  an  equation  of  the  form 


1 


l+ABtj 


+  C 


1 


1  +  CDt 


,  where  A,  B,  C  and  D  are  constants,  to  the  scour 


depth  time  history  from  the  24  hour  scour  experiment  and  then  taking  the  derivative.  The 
form  of  this  equation  is  attractive  because  it  effectively  incorporates  the  two  different 
rates  of  scour  (as  defined  by  Ettema(1980))  and  also  asymptotically  approaches  an 


'  This  form  was  suggested  by  Sterling  Jones  (FHWA),  personal  coiiimunication. 
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equilibrium  value  (A+C).  The  fitted  curve  and  the  data  from  the  experiment  are  plotted 
in  Figure  5.17. 
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Figure  5.17  Experimental  scour  depth  time  history  and  curve  fit  to  the  data. 


Next,  the  KE  data  were  reduced  according  to  the  right  hand  side  of  equation  (5.4). 
Each  record  was  traversed  to  find  the  locations  of  the  excess  kinetic  energy.  This  value 
was  multiplied  by  the  incremental  surface  area  it  covered  and  then  added  to  a  running 
total.  The  best  fit  coefficients  were  then  calculated  using  this  total  and  the  scour  depth 
according  to  equation  (5.4)  in  an  attempt  to  match  it  to  the  curve  produced  using  equation 
(5.2).  The  results  of  the  analysis  are  plotted  in  Figure  5.18. 
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Figure  5.18  Qo  from  equation  (5.2)  and  equation  (5.4)  versus  scour  depth. 


Even  though  there  is  considerable  scatter  in  the  data,  the  results  shown  in  Figure 
5.18  are  encouraging.  The  behavior  of  the  Qo  predicted  using  energy  arguments  resembles 
that  of  the  curve  produced  using  the  experimentally  measured  scour  depth  time  history. 
The  existence  of  the  local  maximum  in  the  curve  of  equation  (5.2)  reflects  two  competing 
processes.  Qo  is  equal  to  a  function  of  order  ds^  multiplied  by  the  time  derivative  of  scour 
depth.  The  rate  of  scour  is  initially  very  high  and  then  rapidly  decreases,  asymptotically 
approaching  zero  while  the  function  of  order  ds^  steadily  increases  from  zero;  hence  the 
existence  of  the  maximum.  Obviously  the  choice  of  a  functional  fit  for  the  scour  depth 
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time  history  plays  a  very  important  role  in  the  form  of  this  curve.  This  is  particularly  true 
in  the  earlier  stages  of  scour  where  most  of  the  scouring  occurs.  The  best  fit  coefficients 
in  equation  (5.4)  for  this  scour  event  are  1.3x10  sec/cm  and  7.6x10^  cm  for  Ci  and  Ci 
respectively.  The  reason  for  the  considerable  scatter  can  be  explained  by  examining  the 
plots  of  the  kinetic  energy  at  the  surface  (Figures  5.12  through  5.15).  Once  the  scour  hole 
has  formed,  FL0W-3D®  predicts  the  existence  of  several  "hot-spots"  on  the  surface  of  the 
holes  where  the  kinetic  energy  experiences  a  local  maximum.  These  "hot-spots"  are 
created  by  the  secondary  flows  inside  the  hole.  The  area  over  which  these  "hot-spots" 
occur  is  highly  dependent  upon  the  scour  depth  and  their  locations  varied  mildly  between 
consecutive  depths. 

To  confirm  the  findings  from  the  simulations,  a  second  set  of  simulations  was  run. 
These  simulations  modeled  a  physical  experiment  that  had  a  different  sediment  size  (dso  = 
0.28  mm),  cylinder  diameter  (b  =  1 1.4  cm),  water  depth  (yo  =  10.4  cm)  and  depth 
averaged  velocity  (U  =  27.3  cm/s).  The  depth  averaged  critical  velocity  was  slightly 
lower  than  in  the  previous  experiments.  A  set  of  5  simulations  was  run  at  depths  equal  to 
2.0, 4.0,  6.0,  8.0  and  9.6  cm.  In  addition,  the  threshold  case  was  run.  Using  the  same 
procedure  as  above,  the  time  history  curve  was  fit  with  a  curve  so  that  the  derivative  may 
be  taken.  The  curve  in  this  case  has  the  same  form  as  before,  namely 


d.  =A 


1   ^  ^ 


,  where  A,B,C  and  D  are  constants.  Figure  5.19 


1+CDt. 

shows  a  plot  of  the  experimental  measurements  and  tlie  best  fit  curve. 


184 


12.0 


10.0 


6.0 


2.0 


- 

  Experimental  Data 

~      d^  =  6.5(l-I/(l  +  1.3xlO"^*t))  +  8.3(l-l/(l  +  1.25xl0"^*t)) 

/ 

1 

/ 

1       1       1       1       1       1       1       1       1       1  1 

O.OE+0       5.0E+3       l.OE+4       1.5E+4       2.0E+4       2.5E+4  3.0E+4 

Time  (seconds) 


Figure  5.19  Scour  depth  time  history  and  best  fit  curve  for  the  second  set  of  simulations. 


The  simulations  were  run  and  the  data  reduced  to  produce  the  plot  in  Figure  5.20. 
The  threshold  simulation  produced  a  KEt  value  of  1 18.0  cm^/s^.  The  lower  threshold 
value  is  a  direct  result  of  the  lower  depth  averaged  critical  velocity.  The  best  fit 
coefficients  for  this  case  are  6.3x10^  sec/cm  and  6.9x10"^  cm  for  Ci  and  C2  respectively. 
The  relatively  large  change  in  the  value  of  Ci  reflects  the  significant  difference  in  the 
values  of  U/Uc,  yo/b,  and  dso/b  for  the  two  tests;  U/Uc  changing  from  0.84  to  0.97;  yo/b 
from  2.1  to  0.91  and  dso/b  from  1.0x10^  to  2.5x10 1  The  relatively  small  change  in  C2 
reflects  the  small  change  in  the  geometry  of  the  hole. 
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Figure  5.20  Qo  from  equation  (5.2)  and  equation  (5.4) 
versus  scour  depth  for  the  second  set  of  simulations. 


10.0 


The  simulation  data  for  this  case  are  not  nearly  as  close  to  the  theoretical  curve  as 
in  the  previous  set.  The  shape  of  the  curve  produced  using  equation  (5.2)  appears  to  be 
highly  dependent  on  the  curve  fit  used  in  the  time  history  data.  Examining  Figure  5.19, 
the  curve  fit  has  a  slightly  different  rate  of  scour  (slope)  at  numerous  points  in  the  record 
than  does  the  experimental  data.  This  shows  up  in  Figure  5.20  as  different  values  of  Qo. 
Equation  (5.2)  proved  to  be  highly  dependent  on  the  form  of  the  equation  chosen  for  the 

^d(ds)^ 


curve  fit  to  the  measured  time  history  data,  since  Qo  is  directly  proportional  to 


dt 


As  a  result,  it  is  important  to  accurately  determine  the  slope  of  this  fit  if  a  meaningful 
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comparison  is  to  be  drawn.  The  time  record  for  this  experiment  is  only  7  hours  long.  A 
longer  experiment  would  most  likely  result  in  a  greater  scour  depth  than  is  predicted  by 
the  fitted  equation.  Larger  scour  depths  would  lead  to  a  greater  slope  in  the  curve  fit  of 
the  time  history  and  thus  a  greater  value  of  Qp  for  the  cofflparison.  In  reality,  the  Qo  for 
this  experiment  is  probably  closer  to  that  predicted  by  the  numerical  simulations  than  that 
developed  from  equation  (5.2).  Regardless,  the  predictions  from  equation  (5.4)  still 
exhibit  the  same  behavior  as  that  of  equation  (5.2).  It  is  expected  that  using  a  longer  time 
history  record  would  improve  the  results  of  the  comparison. 

The  resuUs  from  these  two  sets  of  simulations  indicate  that  this  is  a  promising 
procedure  for  prediction  of  the  depth  of  scour  as  a  function  of  time.  The  area  of  the 
procedure  that  needs  more  work  is  in  determining  the  coefficients  of  proportionality  in 
equation  (5.4).  Arguments  have  been  presented  suggesting  that  the  coefficient  are 
functions  of  the  flow,  structure,  and  sediment  properties  (Ci)  and  the  geometry  of  the  hole 
(C2).  These  hypotheses  are  unconfirmed,  however,  due  to  the  limited  number  of 
simulation  sets  performed  as  part  of  this  investigation.  More  work  is  required  to  find  the 
exact  dependence  of  these  coefficients. 

In  summary,  the  procedure  for  finding  the  depth  of  local,  clear  water  scour  around 
a  circular  cyhnder  as  a  function  of  time  is  as  follows: 

Assumptions:  1)  The  coefficient  of  proportionality  in  equation  (5.4)  is  known  as  a 
function  of  the  flow  and/or  sediment  properties. 

2)  The  scour  hole  can  be  modeled  as  an  inverted  right  cone  with  an  angle 
equal  to  the  angle  of  repose  of  the  sediment.  It  is  truncated  at  the  specific 
depth  of  scour  to  form  a  small  ledge. 
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Procedure:      1)  Using  the  FL0W3D®  software,  run  a  threshold  condition  case  at  a 
depth  averaged  critical  velocity  equal  to  0.45. 

2)  Locate  the  maximum  total  kinetic  energy  on  the  "surface"  at  a 
perpendicular  distance  of  5.0  mm  above  the  bed.  This  is  the  threshold 
kinetic  energy,  KEt,  for  the  subsequent  analysis. 

3)  Using  FL0W3D®,  run  a  number  of  cases  of  varying  depths  from  the 
flat  bed  to  the  scour  depth  where  the  excess  kinetic  energy  is  no  longer 
able  to  exceed  the  effects  of  the  potential  energy  associated  with  the 
removal  of  sediment  from  the  bottom  of  the  hole.  Note  that  this  depth  will 
be  beyond  the  equihbrium  depth. 

4)  From  the  results  of  the  simulations,  integrate  the  excess  kinetic  energy 
in  the  scour  hole  according  to  equation  (5.4). 

5)  Using  die  results  of  4)  and  the  scour  depth  of  the  simulations  to  find  Qo 
as  a  function  of  scour  depth  according  to  equation  (5.4). 

6)  Use  the  results  from  step  5)  in  equation  (5.2)  to  determine  — — ^  as  a 

dt 

function  of  scour  depth,  cylinder  diameter,  ledge  width  and  dynamic  angle 
of  repose  of  the  sediment. 

7)  Numerically  integrate  this  function  in  step  6)  to  determine  scour  depth 
as  a  function  of  time  and  the  equilibrium  scour  depth. 


CHAPTER  6 

SUMMARY,  CONCLUSIONS  AND  RECOMENDATIONS 

6.1  Summary 

With  the  goal  of  providing  a  model  and  methodology  for  predicting  the  temporal 
variation  of  scour,  an  investigation  is  conducted  that  includes  physical  experiments, 
numerical  simulations,  and  the  postulation  of  a  model.  To  better  understand  the  nature  of 
local  scour,  a  thorough  review  of  the  processes  involved  and  the  individual  mechanisms 
present  is  provided.  A  survey  of  the  literature  concerning  both  the  observations  of  the 
rate  of  local  scour  and  current  methods  for  predicting  local  scour  as  a  function  of  time  is 
also  presented.  The  literature  review  concludes  with  a  presentation  of  the  papers  that 
cover  the  measurement  of  hydrodynamics  in  the  scour  hole  made  both  in  the  laboratory 
setting  and  through  the  use  of  computational  experiments. 

Seven  physical  experiments  are  performed  including  three  erodible  bed  scour 
experiments  and  four  fixed  bed  experiments.  The  experiments  were  conducted  in  the 
Hydraulics  Research  Flume  located  in  the  Department  of  Civil  Engineering  at  the 
University  of  Florida.  All  the  experiments  were  conducted  using  a  single  circular  cylinder 
that  is  16.8  cm  in  diameter  placed  in  a  bed  of  sand  that  has  a  median  diameter  of  0. 169 
mm.  The  average  flow  conditions  for  all  the  experiments  are  as  follows:  water  depth  = 
34.8  cm  and  depth  averaged  velocity  =  24.4  cm/s.  The  elapsed  times  for  the 
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erodible  bed  scour  experiments  are  equal  to  one,  six  and  24  hours.  For  each  of  these  tests, 
measurements  were  made  of  the  flow  conditions  and  the  time  history  of  the  depth  of 
scour.  The  four  fixed  bed  hydrodynamic  tests  model  each  of  the  holes  created  during  the 
three  erodible  bed  tests  as  well  as  an  additional  initial  flat  bed  condition.  Measurements  of 
the  velocity  are  made  at  numerous  points  around  the  cylinder  and  in  the  scour  hole  using 
an  Acoustic  Doppler  Velocimeter.  Plots  of  the  velocity  vectors  for  all  of  the 
measurements  are  presented. 

Numerical  simulations  of  each  of  the  fixed  bed  hydrodynamic  tests  were  made 
using  a  three  dimensional,  finite  difference  flow  solver  called  FL0W3D®.  The  sofl;ware  is 
capable  of  handling  complex  geometries  and  can  solve  the  Navier  Stokes  equations  using  a 
variety  of  turbulence  closure  techniques.  The  models  of  the  fixed  bed  tests  use  an 
experimentally  measured  upstream  velocity  profile  as  input.  The  solutions  from  the 
calculations  are  reduced  to  find  the  velocity  vectors  at  the  identical  points  where  the 
measurements  are  made  during  the  physical  experiments.  Plots  of  the  velocity  vectors  at 
these  points  are  presented.  An  additional  set  of  simulations  was  run  using  average 
measured  flow  conditions  from  all  the  fixed  bed  hydrodynamic  tests.  This  set  models  the 
scour  hole  at  20  different  depths  between  the  initial  flat  bed  condition  and  the  scour  hole 
after  24  hours.  A  final  simulation  is  run  using  a  depth  averaged  velocity  equal  to  0.45 
times  the  critical  depth  averaged  velocity  for  the  sediment.  This  condition  is  generally 
agreed  upon  as  the  point  where  scour  is  initiated  for  flow  past  a  single  circular  cylinder. 

The  velocity  vectors  measured  during  the  physical  experiments  are  compared  with 
the  predicted  values  from  the  numerical  simulations.  The  model  accuracy  is  judged  by  the 
ability  to  predict  both  the  measured  velocity  magnitude  and  direction.  Plots  are  presented 
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showing  the  calculated  error  in  the  magnitude  and  direction  as  well  as  difference  vectors 
between  the  model  and  the  experiment.  A  novel  methodology  for  the  prediction  of  the 
temporal  variation  of  scour  is  developed.  Using  the  set  of  21  simulations,  the  excess 
kinetic  energy  at  the  surface  of  the  scour  hole  and  the  potential  energy  associated  with 
removing  sediment  from  the  bottom  of  the  hole  are  found  to  be  proportional  to  the  flux  of 
sediment  leaving  the  scour  hole.  The  results  from  the  simulations  are  presented  with  those 
calculated  theoretically.  To  test  this  method  an  additional  set  of  simulations  modeling  a 
completely  different  scour  situation  is  run.  The  results  from  this  set  are  also  presented. 

6.2  Conclusions 

The  physical  experiments  yielded  a  number  of  interesting  results.  First,  the  method 
for  fixing  the  bed  for  the  hydrodynamic  experiments  was  found  to  be  a  quick,  easy,  yet 
reliable  technique  and  it  is  recommended  that  this  procedure  be  used  for  investigations  of 
this  type  (see  Chapter  3).  The  Sontek  ADV  is  a  sturdy  and  dependable  instrument  that 
provides  repeatable  and  accurate  results.  The  advantages  to  using  this  instrument  includes 
its  ease  of  calibration  and  the  fact  that  it  is  relatively  non-intrusive,  especially  when 
compared  with  such  conventional  methods  such  as  current  balls,  hot-films  or  pitot  tubes. 
The  drawbacks  include  the  fact  that  the  instrument  does  not  perform  well  when  its 
sampling  volume  is  placed  very  close  to  boundaries.  Also,  since  the  cutoff  frequency  is  25 
Hz,  a  significant  portion  of  the  turbulence  for  the  flows  such  as  those  investigated  here  is 
not  being  measured. 

The  measurements  of  the  velocity  vectors  around  the  cylinder  captured  many 
intriguing  features  of  the  flow  field  as  the  hole  is  scoured  out.  In  the  initial  stages  of  scour 
(zero  and  one  hour  scour  holes),  the  main  scouring  mechanisms  are  the  increased  shear 
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stress  due  to  the  acceleration  of  the  flow  around  the  sides  of  the  pile  and  the  downflow  at 
the  leading  edge.  In  the  stagnation  plane,  the  presence  of  a  coherent  horseshoe  vortex  was 
not  observed  for  the  flat  bed  case.  This  may  be  due  to  a  lack  of  sufficient  resolution  of 
measurements  in  the  scour  hole.  For  the  one  hour  scour  hole,  there  was  no  separation  of 
flow  at  the  upstream  edge  of  the  hole.  As  a  result,  the  horseshoe  vortex  measured  was 
quite  small  in  diameter.  It  is  hypothesized  that  the  circulation  associated  with  the 
horseshoe  vortex  does  not  appear  in  the  stagnation  plane  until  there  is  separation  at  the 
upstream  edge  of  the  hole.  In  the  later  stages  of  scour  (6  and  24  hours  scour  holes),  there 
is  separation  at  the  upstream  edge  of  the  hole.  The  horseshoe  vortex  is  easily  identified 
and  it  expands  to  fill  the  entire  cross-section  of  the  hole.  As  it  is  swept  around  the  sides  of 
the  cylinder  it  reduces  in  cross-section  and  it  climbs  up  the  sides  of  the  hole.  There  is  a 
well  defined  circulation  associated  with  the  vortex  even  in  the  90°  plane. 

In  the  wake  region,  the  measurements  show  similar  results  as  the  hole  expands  in 
size.  The  flow  perpendicular  to  the  main  flow  direction  sees  a  marked  increase  in 
magnitude  from  the  flat  bed  condition  to  the  one  hour  scour  hole.  From  this  point,  the 
magnitude  steadily  decreases.  This  may  be  attributed  to  the  fact  that  the  hole  itself 
redirects  fluid  laterally  towards  the  lee  of  the  cylinder.  Upon  further  expansion  of  the 
hole,  the  larger  cross-sections  lead  to  lower  velocities  seen  by  the  pile  and  hence  a  steady 
decrease  in  magnitude.  Examination  of  the  time  records  for  the  wake  region  confirm  the 
results  published  by  Melville  (1975)  that  the  vortex  shedding  decreases  in  frequency  as  the 
scour  hole  deepens. 

The  software  package  FL0W3D®  is  a  well  designed  program  with  a  convenient 
user  interface.  The  ability  to  quickly  design  complex  geometries  proved  to  be  quite  usefiil 
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in  designing  models  of  the  scour  hole  at  various  depths.  Unfortunately,  the  program  in  its 
present,  commercially  available  form  is  not  designed  for  modeling  open  channel  flows.  It 
was  only  after  the  program  was  modified  by  Flow  Science,  Inc.,  that  satisfactory  results 
were  obtained.  The  modification  allowed  the  user  to  input  a  depth  varying  velocity 
(velocity  profile)  at  one  of  the  boundaries.  To  fiarther  improve  the  model's  performance,  a 
method  for  allowing  roughness  at  the  bed,  with  subsequent  turbulence  generation  and 
energy  dissipation,  must  be  developed. 

The  models  of  the  fixed  bed  experiments  showed  the  same  trends  that  were 
measured  during  the  laboratory  experiments.  These  include  the  lack  of  separation  of  flow 
at  the  upstream  edge  of  the  hole  and  the  absence  of  a  significant  horseshoe  vortex  during 
the  initial  phases  of  scour.  In  the  final  two  cases,  separation  occurred  and  a  large 
horseshoe  vortex  that  filled  the  cross- section  of  the  hole  was  found.  Again,  this  vortex 
can  be  seen  narrowing  in  diameter  as  it  is  swept  around  the  sides  of  the  pile  and  it  also 
climbs  up  the  sides  of  the  hole.  In  the  wake  region,  the  same  change  in  the  magnitude  of 
the  lateral  velocities  is  observed,  though  the  differences  are  less  pronounced. 

The  ability  of  the  numerical  simulations  to  predict  the  velocity  vectors  magnitude 
and  direction  was  deemed  capable.  The  model  was  able  to  predict  the  magnitude  of  the 
velocity  vector  within  30%  for  78%  of  the  3906  measurements.  The  three  dimensional 
difference  in  predicted  angle  was  less  than  30°  for  81%  of  the  measurements.  Large 
percentage  errors  in  prediction  of  the  velocity  magnitude  are  attributed  to  the  fact  that  the 
velocities  are  very  low  near  the  boundaries.  Large  differences  in  velocity  direction  are  a 
Sanction  of  the  way  that  the  numerical  model  was  constructed.  The  error  histograms  in 
Chapter  5  indicate  that  there  may  be  a  weak  inverse  relationship  in  the  model's  ability  to 
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accurately  predict  both  magnitude  and  direction.  Although  the  wake  region  intentionally 
lacked  sufficient  care  in  modeling,  the  results  are  not  too  discouraging  in  this  region, 
especially  as  measurements  get  farther  away  from  the  cylinder.  This  was,  however,  the 
area  in  whieh  the  model  performed  the  poorest.  The  symmetry  boundary  condition  for  the 
stagnation  plane  and  the  poor  modeling  of  the  wake  region  contributed  to  weak 
performance  in  this  region. 

Since  the  model  performed  so  well,  it  was  decided  to  use  the  results  from  the 
second  set  of  simulations  in  Chapter  4  to  develop  a  method  for  predicting  the  temporal 
variation  of  scour  as  a  function  of  time.  The  method  is  based  on  energy  arguments 
concerning  the  rate  of  sediment  removal  from  the  scour  hole.  The  flux  of  sediment  leaving 
the  hole  is  hypothesized  to  be  proportional  to  the  difference  in  the  integral  of  the  excess 
kinetic  energy  and  the  potential  energy  required  to  remove  sediment  from  the  bottom  of 
the  hole.  This  relationship  reflects  the  competing  processes  of  the  erosive  and  resistive 
actions  of  the  hydrodynamics  and  the  geometry  of  the  scour  hole.  The  total  kinetic  energy 
was  evaluated  at  a  5  mm  perpendicular  distance  from  the  bed  inside  the  scour  hole.  This 
distance  reflects  the  closest  measurement  to  the  bed  that  can  be  reliably  compared  with  the 
physical  experiments.  The  threshold  condition  (U/Uc  =  0.45)  was  analyzed  to  find  the 
maximum  value  of  the  total  kinetic  energy.  Since  this  is  the  lowest  value  of  depth 
averaged  velocity  ratio  for  which  scour  occurs  around  a  single  circular  cylinder,  the 
maximum  value  of  total  kinetic  energy  found  in  this  case  should  indicate  the  kinetic  energy 
required  to  initiate  sediment  transport. 

The  remaining  20  cases  were  then  examined  to  find  the  areas  in  the  scour  hole 
where  the  total  kinetic  energy  was  greater  than  the  threshold  value.  The  excess  kinetic 
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energy  was  integrated  over  the  areas  where  this  occurs  in  the  scour  hole.  This  quantity 
was  used  along  with  the  scour  depth  to  find  the  best  fit  curve  using  two  free  coefficients. 
The  results  were  found  to  be  proportional  to  the  flux  of  sediment  transport  out  of  the  hole 
when  plotted  against  depth.  The  rate  of  seditiient  transport  out  of  the  hole  was 
determined  by  making  simple  assumptions  about  the  shape  of  the  hole  and  doing  a  curve 
fit  to  the  time  history  of  scour  depth  measured  from  the  experiment.  The  results  were 
encouraging.  To  test  the  hypothesis,  another  set  of  simulations  was  run  using  completely 
different  conditions.  The  results  from  this  set  were  not  as  good  as  the  first  set,  yet 
performed  well  enough  to  warrant  fijrther  work  in  developing  this  technique.  The 
unanswered  question  is  the  nature  of  the  coefficients  of  proportionality  for  the  integrated 
excess  kinetic  energy  and  the  potential  energy.  It  is  hypothesized  that  these  two 
coefficients  are  fianctions  of  the  flow,  structure  and  sediment  parameters  and  the  geometry 
of  the  hole  respectively.  Recommended  improvements  to  this  method  will  be  discussed  in 
the  next  section. 

6.3  Recommendations  for  Future  Work 
The  method  developed  for  determining  the  temporal  variation  of  local  scour  depth 
around  a  single  circular  cylinder  under  clear  water  scour  presents  a  number  of  intriguing 
possibilities.  For  example,  since  the  threshold  kinetic  energy  for  sediment  movement  can 
be  found  by  performing  a  simple  simulation  of  flow  around  a  single  circular  cylinder,  the 
method  should  apply  to  scour  around  more  complicated  objects  such  as  square  piles  or 
pile  groups.  Simulations  of  flow  around  more  complicated  shapes  can  be  easily  performed 
using  the  FLOWS D®  software.  These  simulations  could  be  analyzed  in  the  exact  same 
manner  as  that  used  in  the  circular  cylinder  case.  The  only  additional  requirement  would 
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be  a  reliable  technique  for  modeling  the  shape  of  the  scour  hole  as  it  grows  both  for  use  in 
the  numerical  simulations  and  in  the  calculation  of  the  rate  of  sediment  transport  where 

Qo  =  f(dj|^       J  .  This  technique  could  also  be  used  to  find  the  equilibrium  scour 

depth  for  complicated  shapes  where  empirical  formulae  do  not  exist.  Knowing  the 
threshold  value  from  the  single  circular  cylinder  case,  the  scour  hole  around  a  more 
complicated  shape  could  be  modeled  at  ever  increasing  depths  until  the  predictive  equation 
(equation  (5  .4))  gives  a  zero  value  for  sediment  flux.  This  point  should  equal  the 
equilibrium  depth. 

There  are  a  number  of  drawbacks  to  using  this  technique,  however.  First,  the 
number  of  simulations  required  to  provide  an  adequate  data  set  is  punitive.  At  this  point 
running  a  single  simulation  for  the  simple  cases  studied  here  required  between  12  to  24 
hours.  A  set  of  20  simulations  would  require  almost  a  month  of  work.  However,  with 
every  increase  in  computer  speed  and  with  every  improvement  in  the  software's  code,  this 
time  will  be  reduced.  Also,  until  the  software  is  modified  to  allow  both  depth  varying 
velocity  input  and  the  problem  of  roughness  at  the  boundaries  is  resolved,  this  will 
continue  to  be  a  deficiency  in  the  simulations. 

Future  work  on  this  method  should  concentrate  on  obtaining  a  relationship 
between  the  energies  associated  with  the  erosion  of  sediment  and  the  rate  of  sediment 
transport.  The  coeflficients  of  proportionality  appear  to  be  a  fianction  of  the  flow 
properties  and  the  geometry  of  the  hole,  however  this  is  unconfirmed.  Investigations 
should  be  pursued  to  find  the  true  nature  of  these  coeflficients. 
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In  addition,  the  inclusion  of  the  stochastic  nature  of  this  process  should  be 
considered.  Since  both  the  hydrodynamics  and  sediment  transport  are  stochastic 
processes,  the  probabilities  assigned  to  different  scour  depths  should  be  investigated.  This 
type  of  information  would  be  quite  usefiil  in  bridge  design  and  bridge  monitoring. 

Also,  improvements  can  be  made  in  the  calculation  of  the  threshold  kinetic  energy. 
It  should  be  related  to  the  critical  depth  averaged  velocity  or  at  least  the  critical  shear 
stress.  An  empirical  relation  for  determining  this  value  would  reduce  by  one  the  number 
of  simulations  required  for  the  implementation  of  this  method. 
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Program  1-  "dataredu.m" 
"/oDATAREDU.M  -  This  infile  will  reduce  the  data  taken  using  the  sontek  acoustic  velocimeter. 
%Many  of  the  variables  required  for  this  reduction  are  listed  in  the  file  sumvar.dat  that 
%is  contained  in  each  of  the  subdirectories  listed  below.  For  a  listing  of  the  organization 
%of  this  file  see  the  text  document  TEMPLATE.TXT  listed  at  the  end  of  this  program  listing. 
%The  mfile  GETANG.M  is  used  in  this  program  and  is  also  listed  at  the  end  of  this  program. 

clear  variables; 

path(path,'c:\sontek');  %add  getang  to  path 

k=fflenu('Choose  a  case:',TlatbedVScouratt  =  tl','Scouratt  =  t2','Equilibrium  Scour'); 
ifk==l, 

eval(['cd  flatbed']); 

end 

if  k  ==  2, 

eval(['cdtl']); 

end 

ifk==  3, 

eval(['cd  t2']); 

end 

ifk==  4, 

eval(['cd  ultscour']); 

end 

fidO=fopen('sumvar .  dat') ; 


p=str2num(fgeU(fidO)); 

xzvO=str2num(fgetl(fidO)); 

yzvO=str2num(fgetl(fidO)); 

shiftilt(l)-str2num{fgetl(fidO)); 

shiftiU(2)=str2num(fgetl(fidO)); 

shiftiU(3)=str2num(fgetl(fidO)); 

shiftref(  1  )=str2num(fgefl(fidO)); 

shiftref(2)=str2num(fgetl(fid0)); 

shiftref(3)=str2num(fgetl(fidO)); 


%the  number  of  test  subdirectories  (days) 
%base  voltage  in  vertical  position  (green  wire) 
%base  voltage  in  vertical  position  (white  wire) 
%X  position  change  vetor  due  to  probe  tilt 
%Y  posiUon  change  vetor  due  to  probe  tilt 
%Z  position  change  vetor  due  to  probe  tilt 
%X  zero  reference  to  center  of  cylinder 
%Y  zero  reference  to  center  of  cylinder 
%Z  zero  reference  to  center  of  cylinder 


ftot  =  0  %total 
for  dirnum  =  1  :p, 

linefile=fgetl(fidO); 

day  =  linefile(l:2); 

fl  =  1; 

f2  =  str2num(linefile(4:6)); 
q  =  str2num(linefile(8:10)); 
flt=str2num(linefile(12: 14)); 
f2t=str2num(linefile(16: 18)); 
xzvol=str2num(linefile(20 : 24)) ; 
yzvol=str2num(linefile(26 : 30)); 
eval(['cd '  day]); 


for  given  test  configuration  (running  variable) 


%test  subdirectory  by  day 

%the  first  file  number  in  the  subdirectory 

%the  final  file  number  in  the  subdirectory 

%missing  file  number 

%start  file  number  for  tilted  position 

%end  file  number  for  tilted  position 

%green  wire  voltage  for  tilted  position 

%white  wire  voltage  for  tilted  position 


fori  =  fl:f2; 

ftot  =  ftot  +  1; 
ifi==  q, 
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outfile(ftot,l)  =  dimum; 
%outfiIe  will  contain  the  summary  of  the  data  from  all  the  reduction 
outfile(ftot,2)  =  i 

elseif  i  ~=  q,  %filename  contains  the  name  of  the  active  file 

ifi<  10, 

filename  =  [  day  '00'  num2str(i)  ]; 
elseif  i  >=  10, 

ifi<  100, 

filename  =  [  day  '0'  num2str(i)  ]; 
elseif  i>=  100, 

filename  =  [  day  num2str(i)  ]; 

end 

end 

eval(['!..\..\getctl '  filename]); 
eval{['!..\..\getvel '  filename]); 
eval(['load '  filename  '.vel']); 
infile=eval(filename) ; 
eval(['!del '  filename  '.vel']); 
eval(['clear '  filename]); 


%ctl  file  manipulation 


ctlfile= [filename  '.ctl']; 
fid=fopen(ctlfile); 
for  j= 1:4, 

dum=fgetl(fid); 

end 

linex=fgetl(fid); 

liney=fgetl(fid); 

linez=fgetl(fid); 

nx=length(linex); 

ny=length(liney); 

nz=length(linez); 

x=str2num(linex(  1 6 :  nx)) ; 

y=str2num(liney(  IS.ny)); 

z=str2num(linez(18:nz)); 

forj=l;19, 

dum=fgetl(fid); 

end 

linezbot=fgetl(fid) ; 
iflinezbot(ll)=^  'N', 

zpt=999; 

zsv=999; 


%ctlfile  is  the  name  of  the  .ctl  file 
%that  is  generated  by  the  adf  software 


%x  position  from  sontek  input 
%y  positiOit  from  sont6k  input 
%z  position  from  sontek  input 


else 


%zpt  is  distance  from  probe  tip  to  bottom 
%zsv  is  distance  from  sample  volume  to  bottom 


zpt=str2num(linezbot(  15:19)); 

%zpt  and  zsv  are  given  a  value  of  999  if 
zsv=str2num(linezbot(36:40)); 

%the  bottom  is  'N'ot  found. 


end 

fclose(fid); 
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eval(['!del '  filename  '.ctl']); 


m=length(infile( : ,  1 )) ; 

pos=[x  y  z];  %  (x,y,z)  position  vector 

% 

%Probe  rotation  Section 

% 

if  xzvol-^  9995^9, 

xztO=-getang(  1  ,xzvol,xzvO);  %xztO  and  yztO  are  the  angles  that  the  probe  is 
yzt0=-getang(2,yzvol,yzv0);     %rotated  as  determined  by  the  mfile  getang.m 

xztO=xztO*2*pi/360; 

yztO=yztO*2*pi/360; 

Rx  =  [1  0  0;  0  cos(yztO)  sin(yztO);  0  -sin(yztO)  cos(yztO)  ]; 

%Rx  =  rotation  matrix  in  xz  plane 
Ry  =  [cos(xztO)  0  -sin(xztO);  0  1  0;  sin(xztO)  0  cos(xztO)  ]; 

%Ry  =  rotation  matrix  in  yz  plane 

A  =  Rx*Ry; 

%Rotation  matrix  from  tilted  to  absolute  coordinates 

% 

%Reference  for  these  calculations  is  Mathematical  Methods  for  Physicists  by  Arfken  1970,  pgs  178-79 
% 

ifi>=  fit, 

ifi< 


end 

end 

end 

% 

%Calculation  of  output  section 

0/ 

/o 

pos  =  pos  +  shiftref;  ^Reference  position  vector  to  center  of  cylinder 
pos(3)  =  -(pos(3));  %and  convert  all  units  to  cgs 

pos(l)  =  pos(I)*100; 
pos(3)  =  pos(3)/10; 


=  f2t, 

for  j  =  l:m, 

Vt=[infile(j,4)  infile(j,5)  infile(j,6)]; 
%Velocity  matrix  in  tilted  coordinates 
V  =  A*Vt'; 

%Velocity  matrix  corrected  for  rotation 

infile(j,4)=V(l); 

infile(j,5)=V(2); 

infile(j,6)=V(3); 

end 

pos  =  pos  +  shiftilt; 

"/ocorrect  position  for  tilt 


u=infile(:,4); 
v=infile(:,5); 
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w=infile(:,6); 

% 

%  In  this  section,  the  data  is  fiUered  first  to  eliminate  the  single  point  spikes  that 

%  tend  to  occur  in  sontek  data  and  secondly  to  remove  the  low  frequency  (below  0.2  Hz) 

%  oscillation  that  is  found  in  the  flume.  The  spikes  are  removed  by  finding  all  excursions 

%  that  are  more  than  3  standard  deviations  from  the  previous  data  point  and  setting  them  equal 

%  to  the  mean  of  the  bracketing  points.  The  low  frequency  oscillation  is  removed  by  taking  the 

%  fourier  transform  of  the  time  series  and  setting  equal  to  zero  all  the  points  that  are 

%  below  the  0.2  Hz  cutoff.  Then  the  iffl  is  taken  to  recover  the  time  series. 

% 

su=std(u); 
sv=std(v); 
sw=std(w); 

ut(l)=u(l); 

ut(m)=u(m); 

vt(l)=v(l); 

vt(m)=v(m); 

wt(l)=w(l); 

wt(m)=w(m); 

fork=2:m-l, 

ut(k)=u(k); 

dul=abs(u(k)-ut(k-l)); 
du2=abs(u(k+l)-ut(k-l)); 
if  dul  >  3*su, 

ifdul>2*du2, 

ut(k)=0,5*(u(k+l)+ut(k-l)); 

end 

end 

vt(k)=v(k); 

dvl=abs(v(k)-vt(k-l)); 
dv2=abs(v(k+l)-vt(k-l)); 
if  dvl  >  3*sv, 

ifdvl  >2*dv2, 

vt(k)=0. 5  *(v(k+l  )+vt(k-l )); 

end 

end 

wt(k)=w(k); 

dw  1 =abs(w(k)-wt(k- 1 )); 
dw2=abs(w(k+ 1 )  -wt(k- 1 )) ; 
if  dwl  >  3*sw, 

if  dwl  >  2*dw2, 

wt(k)=0.5*(w(k+l)+wt(k-l)); 

end 

end 

end 

a=ceil(0.2*infile(m,l)); 

b-m-a; 

fii=fft(ut); 

fu(2  :a)=zeros(size(2  :a)); 
fu(b  :m)=zeros(si  ze(b :  m)) ; 
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u=real(iffl(fu)); 
fV=ffl(vt); 

fV(2 :  a)=zeros(size(2 :  a)) ; 
fV(b :  m)=zeros(size(b :  m)); 
v=real(iffl(fV)); 
fw=ffi(wt); 

fw(2 :  a)=zeros(size(2 :  a)); 
fw(b :  m)=zeros(size(b :  m)) ; 
w=real(ifit(fw)); 

% 

%  Calculation  of  flow  statistics 

% 

um  =  mean(u); 

vm  =  mean(v); 

wm  =  mean(w); 

up=u  -  um; 

vp=v  -  vm; 

\vp=w  -  wm; 

upupm  =  mean(up.*up); 

vpvpm  =  mean(vp.*vp); 

wpwpm  =  mean(wp.*wp); 

upvpm  =  mean(up.*vp); 

upwpm  =  mean(up.*wp); 

vpwpm  =  mean(vp.*wp); 

s  =  [upupm  upvpm  upwpm;  upvpm  vpvpm  vpwpm;  upwpm  vpwpm  wpwpm]; 
d  =  det(s); 

outfile(ftot,l)  =  dimum; 
outfile(ftot,2)  =  i; 
[dimum,i] 

outfile(ftot,3)  =  infile(m,l); 
outfile(ftot,4)  =  um; 
outfile(ftot,5)  =  vm; 
outfile(ftot,6)  =  wm; 
outfile(ftot,7)  =  pos(l); 
outfile(ftot,8)  =  pos(2); 
outfile(ftot,9)  =  pos(3); 
outfile(ftot,10)  =  zpt; 
outfile(ftot,l  1)  =  zsv; 
outfile(ftot,12)  =  upupm; 
outfile(ftot,  13)  =  vpvpm; 
outfile(ftot,14)  =  wpwpm; 
outfile(ftot,15)  =  upvpm; 
outfile(ftot,16)  =  upwpm; 
outfile(ftot,17)  =  vpwpm; 
outfile(ftot,18)  =  d; 

end 

end 

eval(['cd  ..']); 

end 

save  totalout.dat  outfile  -ascii 

fclose(fidO) 

eval(['cd ..']); 
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function  theta=getang(xy,vin,vO) 

%GETANG.M  This  Mfile  calculates  the  angle  that  the  inclinometer  is  tilted 
%from  the  given  voltages.  The  fiinctions  are  a  best  fit  5th  order  polynomial 
%that  was  supplied  by  the  manufacturer:  APPLIED  GEOMECHANICS  INC. 


a(l)=0; 

b(l)=41.7924; 

c(l)=-0.6146; 

d(l)=-5.7194; 

e(l)=0.0329; 

f(l)=0.6567; 

a(2)=0; 

b(2)=4 1.7886; 

c(2)=-1.5157; 

d(2)=-5.7329; 

e(2)=0,2582; 

f(2)=0.6636; 


thetO=a(xy)  +  b(xy)*vO  +  c(xy)*v0'^2  +  d(xy)*vO'^3  +  e(xy)*v0'^4  +  f(xy)*vO'^5; 
thetirl=a(xy)  +b(xy)*vin  +  c(xy)*vin'^2  +  d(xy)*vin'^3  +  e(xy)*vin'^4  +  f(xy)*vin'^5; 
theta  =  thetin  -  thetO; 


TEMPLATE.TXT:  Data  found  in  sumvar.dat  has  the  following  organization: 


p  the  number  of  test  subdirectories  (days) 

xzvO  base  voltage  in  vertical  position  (green  wire) 

yzvO  base  voltage  in  vertical  position  (white  wire) 

shiftilt(  1 )  X  position  change  vetor  due  to  probe  tilt 

shiftilt(2)  Y  position  change  vetor  due  to  probe  tilt 

shiftilt(3)  Z  position  change  vetor  due  to  probe  tilt 

shiftref(l)  X  zero  reference  to  center  of  cylinder 

shiftref(2)  Y  zero  reference  to  center  of  cylinder 

shiftref(3)  Z  zero  reference  to  center  of  cylinder 

day       f2  q  fit        fZt        xzvol  yzvol 


where 

day  =  test  subdirectory  by  day  (e.g.  fO,  bl  etc.) 

f2  =  the  final  file  number  in  the  subdirectory 

q  =  missing  file  number 

fit  =  start  file  number  for  tilted  position 

f2t  =  end  file  number  for  tilted  position 

xzvol  =  green  wire  voltage  for  tilted  position 

yzvol  =  white  wire  voltage  for  tilted  position 


Note,  if  xzvol  and  yzvol  =  99999  then  probe  is  in  vertical  position  for  all  runs  on  this  day. 
File  naming  convention: 
e.g.  al095.adv 

First  two  characters  (al)  correspond  to  the  day  (for  specific  days  of  tests  reference  lab 

notebook). 
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Program  2  -  "compare. for" 

*  * 

*  Program  compare.for  * 

*  * 

*  This  program  reads  both  the  output  files  created  by  the  M-file  * 

*  dataredu.m  and  the  Ascii  output  from  the  software  FL0W3D.  ft  * 

*  then  uses  a  three  dimensional  interpolation  scheme  to  determine  * 

*  the  velocity  predicted  by  the  flow  solver  at  the  same  point  * 

*  where  the  measurement  was  made  during  the  physical  experiments.  * 

*  The  interpolation  scheme  is  a  subroutine  located  in  the  IMSL  * 

*  libraries.  * 

*  + 

*  Input  Files:        flowinfo.dat       (fromFLOW3D)  * 

*  totout.dat  (from  datiircdu.m)  * 

*  Output  Files:      compare.dat  * 

*  * 

*  Written  by  Mark  Gosselin  (1/10/96)  * 

*  * 

*  Last  Modified:  5/20/97  * 

*  + 


program  compare 
use  msimslms 

real  xc,yc,zc,uc,vc,wc,fc,xchalf,ychalf,zchalf,sontek 
real  xs,ys,zs,us,vs,ws,uinterp,vinterp,winterp 
real  velcmag,velsmag,vmagerr,vdirerr,vdot 
integer  icjc,kc,nic,njc,nkc,is,nis,js 
character*98  dumchar 

dimension  xc(60),yc(60),zc(60),uc(60,60,60),vc(60,60,60) 
dimension  wc(60,60,60),fc(60,60,60),xchalf(60),ychalf(60) 
dimension  zchalf(60),sontek(  1 500, 1 8) 
logical  check 


*  Data  Input  Section 

* 

open(l,file='flowinfo.dat') 
rewind(I) 

read(l,*)  dumchar 
10        if  (.not.eof(l))then 

read(l,*)icjc,kc,xc(ic),yc(jc),zc(kc),uc(icjc,kc), 
&  vc(ic,jc,kc),wc(icjc,kc),fc(ic,jc,kc) 
print  *,ic,jc,kc 
goto  10 

endif 
close(l) 


nic=ic 
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njc=jc 
nkc=kc 

+ 

*  Note-the  flow  velocities  determined  by  FL0W3D  are  located  at  the 

*  center  of  the  corresponding  cell  faces. 
* 

do  20  ic=2,nic 

xchalf(ic- 1  )=(xc(ic)+xc(ic- 1  ))/2 . 

20  continue 

do  30  jc=2,njc 

ychalfljc- 1  )=(yc(jc)+yc(jc- 1  ))/2 . 

30  continue 

do  40  kc=2,nkc 

zchalf(kc- 1  )=(zc(kc)+zc(kc- 1  ))/2 . 

40  continue 

open(2,file='totalout.dat') 
rewind(2) 

is=0 

50        if(.not.eof(2))  then 

read(2,*)  (sontek(is,js),js=l,18) 

is=is+l 

goto  50 

endif 
close(2) 

nis  =  is  -  1 
check  =  .true. 

*  Interpolation  of  Data 
* 

open(3,file='compare.dat',status='unknown') 

rewind(3) 

do  100  is  =  l,nis 

xs  =  sontek(is,7) 

ys  =  sontek(is,8) 

zs  =  sontek(is,9) 

us  ^  sontek(is,4) 

vs  =  sontek(is,5) 

ws  =  sontek(is,6) 

xs  ~  -xs 

ys  =  -ys 

us  =  -us 

vs  =  -vs 
print  *,'is  =  ',is 

if((xs.eq.0.0).and.(ys.eq.0.0).and.(zs.eq.0.0).and. 
&         (us.eq.0.0).and.(vs.eq.O.0),and.(ws.eq.0.O))  then 
uinterp  =  0.0 
vinterp  =  0.0 
vvinterp  =0.0 
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vmagerr  =  0.0 
vdirerr  =  0.0 
goto  90 

endif 

uinterp  =  QD3VL(xs,ys,zs,mc,xc,njc-l,ychalf,nkc-l,zchalf, 
&  uc(  1  :nic,2  :njc,2  :nkc),nic,njc- 1  ,check) 

vinterp  =  QD3VL(xs,ys,zs,mc-l,xchalf,njc,yc,nkc-l,zchalf, 
&  vc(2:nic,  1  :njc,2:nkc),nic-l,njc,check) 

winterp  =  QD3VL(xs,ys,zs,nic-l,xchalf,njc-l,ychalf,nkc,zc, 
&  wc(2  :nic,2:njt,  1  :iikc),nic- 1  ,njG- 1  ,check) 

velcmag  =  (uinterp**2.  +  vinterp**2.  +  winterp**2.)**0.5 

velsmag  =  (us**2.  +  vs**2.  +  ws**2.)**0.5 

vmagerr  =  abs((velcmag  -  velsmag)/velsmag)*100. 

vdot=us*uinterp  +  vs*vinterp  +  ws*wiiiterp 

if  ((velcmag*velsmag).lt.0.01)  then 
vdirerr  =  999.0 

endif 

if  ((velcmag*velsmag).ge.0.01)  then 

vdirerr  =  acos(vdot/(velcmag*velsmag)) 

endif 

90  write(3,1004)  xs,ys,zs,us,vs,ws,uinterp,vinterp,winterp, 

&  vmagerr,vdirerr 
100  continue 
* 

*  Note  -  the  error  between  the  measured  and  predicted  velocities 

*  is  represented  by  the  error  in  magnitude  (vmagerr)  and  in 

*  direction  (vdirerr). 


1000  format(a98) 

1001  format(3i2,7e7.5) 

1002  format(llal4) 
1004  format(llel4.6) 

end 
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Program  3  -  "energy.for" 

* 

*  Program  energy.for 
* 

*  This  program  uses  the  output  ASCII  files  fi-om  the  FLOW3D 

*  software.  It  uses  a  three  dimensional  interpolation  scheme  to 

*  calculate  the  velocities  and  turbulent  kinetic  energy  at  points 

*  located  at  a  5  mm  perpendicular  distance  fi-om  the  surface  of  the 

*  scour  hole.  The  calculation  mesh  is  80x40  cm  in  the  x  and  y 

*  direction  respectively  at  a  spacing  of  5  mm.  The  program  makes 

*  use  of  the  IMSL  subroutine  QD3VL  for  the  3d  interpolation,  ft  also 

*  uses  the  subroutine  fmdz  located  at  the  end  of  this  program. 


* 

*  Input  files:  turbinfo.dat 

*  flowinfo.dat 

*  ds.dat  * 

*  ♦ 

*  Output  file:        energy.dat  * 

*  * 

*  Subroutines:       fmdz  * 

*  QD3VL  (IMSL)  * 

*  * 

*  Created  by:        Mark  Gosselin  (2/19/97)  * 

*  Last  Modified:     5/24/97  * 


program  energy 
use  msimslms 

realxc,yc,zc,uc,vc,wc,fc,xchaU;ychalf,zchalf,q,dtke 
real  xs,ys,zs,uinterp,vinterp,winterp,ds 
real  ke,xt,yt,qinterp 
integer  ic,jc,kc,nic,njc,nkc,ixt,iyt 
character*98  dumchar 

dimension  xc(60),yc(60),zc(60),uc(60,60,60),vc(60,60,60) 
dimension  wc(60,60,60),fc(60,60,60),xchalf(60),ychalf(60) 
dimension  zchalf(60),q(60,60,60),dtke(60,60,60) 
logical  checkjncyl 


Data  Input  Section 

open(4,file='ds.dat') 
rewind(4) 
read(4,*)  ds 
close(4) 

open(2,file='turbinfo.dat') 
rewind(2) 

print  *, 'Processing  turbinfo.dat' 
read(2,*)  dumchar 
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5  if(.not.eof(2))then 

read(2,*)icjc,kc,xc(ic),yc(jc),zc(kc),q(ic,jc,kc), 
&  dtke(icjc,kc) 
goto  5 

endif 
close(2) 


open(  1  ,file='flowinfo.dat') 
rewind(i) 

print  *,'Processing  flowinfo.dat' 
read(l,*)  dumchar 
10  if(.not.eof(l))then 

read(l,*)icjc,kc,xc(ic),yc(jc),zc(kc),uc(ic,jc,kc), 
&  vc(ic,jc,kc),wc(ic,jc,kc),fc(icjc,kc) 
goto  10 

endif 

close(l) 

nic=ic 
njc=jc 
nkc=kc 

* 

*  Note  that  FL0W3D  evaluates  the  velocities  at  the  faces  of  the  control  volumes 

*  and  the  turbulent  kinetic  energy  at  the  center  of  the  control  volume. 
* 

do  20  ic=2.nic 

xclialf(ic- 1  )=(xc(ic)+xc(ic- 1  ))/2 . 

20  continue 

do  30  jc=2,njc 

ychalf(jc-  l)=(yc(jc)+yc(jc- 1  ))/2 . 

30  continue 

do  40  kc=2,nkc 

zchalf(kc- 1  )=(zc(kc)+zc(kc- 1))/2. 

40  continue 

check  =  .true. 

*  Data  processing  section  -  output  is  contained  in  the  file  energy.dat 
* 

open(3,file='energy.dat',status='unknown') 
rewind(3) 


do  200  ixt  =  -80,80 

xt  =  float(ixt)/2.0 
do  100  iyt  =  -80,0 

yt  =  float(iyt)/2.0 

incyl=. false. 

call  findz(xt,yt,ds,incyl,xs,ys,zs) 
if  (incyl)  goto  100 

uinterp  =  QD3VL(xs,ys,zs,nic,xc,njc-l,ychalf,nkc-l,zchalf, 
&  uc(I  :nic,2:njc,2:nkc),nic,njc-],check) 
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vinterp  =  QD3VL(xs,ys,zs,nic-l,xchalf,njcjc,nkc-l,zchalf, 
&  vc(2:nic,l:njc,2:nkc),nic-l,njc,check) 

winterp  =  QD3VL(xs,ys,zs,nic-l,xchalf,njc-l,ychalf,nkc,zc, 
&  wc(2  :nic,2  :njc,  1  :nkc),nic- 1  ,njc- 1  ,check) 

qinterp  =  QD3VL(xs,ys,zs,nie-l,xehalCiijG-l,ychalf,nkc-l, 
&  zchalf,q(2  :nic,2 :  njc,2  :nkc),nic- 1  ,n)c- 1 , check) 

ke  =  uinterp**2.0  +  vinterp**2.0  +  winterp**2.0 
ke  =  ke*0.5  +  qinterp 
90  write(3,1004)  -xt,-yt,ke,-xs,-ys,zs 

100  continue 

print  *,-xt 
200  continue 


1000  fonnat(a98) 

1001  fomiat(3i2,7e7.5) 

1002  format(llal4) 
1004  format(llel4,6) 

end 


subroutine  findz(xt,yt,ds,incyl,xs,ys,zs) 

This  subroutine  finds  the  point  5  mm  perpendicular  from  the  surface  of  the 
scour  hole  for  the  given  input  location  (xs,  ys)  on  the  calculation  mesh. 

real  xt,yt,zt,ds,xs,ys,zs,phi,rc,rl,r,theta,pi,dx,dz,dy4elta 
logical  incyl 


rc  =  0.5*6.625*2.54 

rl=  1.0 

delta  =  0.5 

pi  =  4.0*atan(1.0) 

phi  =  31.22*pi/180 

r  =  sqrt(xt**2.0  +  yt**2.0) 

if  (r.le.rc)  then 

incyl  =  true. 

return 

endif 


if  (ds.eq.0.0)  then 
xs  =  xt 
ys  =  yt 
zs  =  delta 
return 

endif 


if  (ds.ne.0.0)  then 

if  (r.h.(rc+rl))then 
xs  =  xt 
ys  =  yt 

zs  =  delta  -  ds 
return 

endif 


210 


if  (r.gt.(rc+rl+ds/tan(phi)))  then 
xs  =  xt 
ys  =  yt 
zs  =  delta 
return 

endif 

zt  =  r*tan(phi)  -  (ds  +  tan(phi)*(rc  +  rl)) 
dz  =  delta*cos(phi) 

if(xt.eq.O.O)then 
dx  =  0.0 

dy  =  delta*sin(phi) 

endif 

if  (xt.ne.O.O)then 

theta  =  atan(yt/xt) 

dx  =  delta*sin(phi)*abs(cos(theta))*(xt^abs(xt)) 
dy  =  delta*sin(phi)*abs(sin(theta)) 

endif 

zs  =  zt  +  dz 
xs  =  xt  -  dx 
ys  =  yt  -  dy 
return 

endif 

return 

end 
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